Rosenbrock and linearly implicit methods

Let $D\subseteq \mathbb{R}\times {\mathbb{R}}^d$ be open, let $d\ge 1$, and let $f:D\to {\mathbb{R}}^d$ be a vector field (First-order systems, initial value problems, and solutions on intervals). For $(t,y)\in D$, a step size $h>0$, a parameter $\gamma >0$, and a matrix $\tilde{J}\in {\mathbb{R}}^{d\times d}$, put $A_{\gamma ,h}(\tilde{J}):=I_d-\gamma h\tilde{J}\mathrm{.}$ When this matrix is invertible, the one-stage Rosenbrock step uses the stage slope $k$ and update $A_{\gamma ,h}(\tilde{J})k=f(t,y),{\mathrm{\Phi }}_{\gamma ,h}^{\tilde{J}}(t,y):=y+hk\mathrm{.}$ At points where all state partial derivatives exist, the exact-Jacobian branch sets $\tilde{J}=D_{y}f(t,y)$, where $D_{y}f$ is the Jacobian in the state variables as in The Jacobian matrix of partial derivatives and the gradient in the scalar-valued case. An approximate-Jacobian branch may supply a different $\tilde{J}$; the equation requires the chosen shifted matrix to be invertible. The name Rosenbrock–Euler here refers to the case $\gamma =1$. This definition specifies one linearized step and does not say that it solves the full backward-Euler nonlinear equation.

Let $d\ge 1$, $B\in {\mathbb{R}}^{d\times d}$, $h>0$, $\gamma >0$, and $0\le q<1$ satisfy $\mathrm{\parallel }\gamma hB{\mathrm{\parallel }}_2\le q$, where $\mathrm{\parallel }\cdot {\mathrm{\parallel }}_2$ is the matrix norm induced by the Euclidean vector $2$-norm, as in The matrix norm induced by a published vector p-norm. Then $A=I_d-\gamma hB$ is invertible and $\mathrm{\parallel }{A}^{-1}{\mathrm{\parallel }}_2\le \frac{1}{1-q}\mathrm{.}$ In particular, if $\mathrm{\parallel }\gamma hB{\mathrm{\parallel }}_2\le \frac{1}{2}$, then $\mathrm{\parallel }{A}^{-1}{\mathrm{\parallel }}_2\le 2$.

Let $d\ge 1$, $h,\gamma >0$, and $J,\tilde{J}\in {\mathbb{R}}^{d\times d}$. Set $A=I_d-\gamma hJ$ and $\tilde{A}=I_d-\gamma h\tilde{J}$. If both matrices are invertible, then ${\tilde{A}}^{-1}-A^{-1}=\gamma h{\tilde{A}}^{-1}(\tilde{J}-J)A^{-1}\mathrm{.}$ Consequently, if $\mathrm{\parallel }{A}^{-1}{\mathrm{\parallel }}_2\le K$ and $\mathrm{\parallel }{\tilde{A}}^{-1}{\mathrm{\parallel }}_2\le \tilde{K}$, then $\mathrm{\parallel }{\tilde{A}}^{-1}-A^{-1}{\mathrm{\parallel }}_2\le \gamma hK\tilde{K}\mathrm{\parallel }\tilde{J}-J{\mathrm{\parallel }}_2\mathrm{.}$ Here the induced matrix norm is the one defined in A small shifted Jacobian has an explicit inverse bound.

Let $d\ge 1$, $M\in {\mathbb{R}}^{d\times d}$, and $x,b\in {\mathbb{R}}^d$ be stored real data. Use successful sequential left-to-right dot products and componentwise subtraction in the operationwise model of Operationwise real and complex floating point arithmetic with rounded square roots and fused multiply-add, with unit roundoff $u>0$ and $du<1$. Put ${\widehat{v}}_i:=\mathrm{fl}\left({\sum }_{j<d}{M}_{ij}{x}_j\right),\widehat{r}:=\mathrm{fl}(\widehat{v}-b),{\gamma }_d:=\frac{du}{1-du}\mathrm{.}$ For $\mathrm{\parallel }M{\mathrm{\parallel }}_F=({\sum }_{i,j<d}{M}_{ij}^2{)}^{1\mathrm{/}2}$, define $W:={\gamma }_d\mathrm{\parallel }M{\mathrm{\parallel }}_F\mathrm{\parallel }x{\mathrm{\parallel }}_2+u(\mathrm{\parallel }\widehat{v}{\mathrm{\parallel }}_2+\mathrm{\parallel }b{\mathrm{\parallel }}_2)\mathrm{.}$ Then $\mathrm{\parallel }\widehat{r}-(Mx-b){\mathrm{\parallel }}_2\le W,\mathrm{\parallel }Mx-b{\mathrm{\parallel }}_2\le \mathrm{\parallel }\widehat{r}{\mathrm{\parallel }}_2+W\mathrm{.}$ If a norm in these formulas is replaced by an exact or certified upper bound, the resulting value of $\mathrm{\parallel }\widehat{r}{\mathrm{\parallel }}_2+W$ remains a certified upper bound for the exact residual.

Input: a stored real square matrix $M\in {\mathbb{R}}^{d\times d}$ with $d\ge 1$; a unit roundoff $u>0$ with $du<1$; and a terminating certificate evaluator evaluator. For each requested Euclidean or Frobenius norm of stored data, evaluator returns a finite stored nonnegative upper bound for the exact norm or unavailable. For each requested nonnegative expression formed from such upper bounds, $u$, and $d$, it returns a finite stored nonnegative upper bound for the exact expression or unavailable. Every returned bound is sound. All reflector, substitution, dot-product, and residual operations below use Operationwise real and complex floating point arithmetic with rounded square roots and fused multiply-add; certificate-evaluator calls return certified bounds rather than relying on ordinary rounded evaluation of a bound expression.

Denote the top-level factor-only operation by factor(M,u,evaluator). It returns either a reusable factor object or a factorization failure status.

1. Set $R:=M$. For each column $j=0,\dots ,d-1$, let $x$ be the active part of column $j$ of $R$ and $B$ the active trailing block to its right. Call Scaled Householder reflector formation and application on $(x,B)$. If it returns range-failure, return that status and no factor object. If it returns a marked identity or reflector result, replace the active block by its returned $C_{+}$, store its marker and (when present) $(v,\widehat{\tau })$, and continue. For any other return, stop and return factorization-failure with no factor object. The resulting $R$ is upper triangular with its lower entries set to exact zeros by the called procedure; this is the full/reduced computational QR convention of Full, reduced, and column-pivoted computational QR factorisations.
2. If any diagonal entry $R_{jj}$ is zero, return singular-factor. Otherwise return a factor object containing the original stored $M$, $R$, the ordered stored reflectors, $u$, and the certificate evaluator. The factor object's operation solve(b,epsilon), whose input is any stored real $b\in {\mathbb{R}}^d$ and finite $ϵ\ge 0$, is:
   1. Apply the stored reflector updates in the same order, each using the successful rounded dot, scaling, product, and subtraction operations specified by Scaled Householder reflector formation and application, to obtain $c$. Return the first range-failure if an operation fails.
   2. Compute $x$ by descending rounded back substitution for $Rx=c$. For $i=d-1,d-2,\dots ,0$, initialize $s:=0$; for $j=i+1,i+2,\dots ,d-1$ in increasing order, compute $p:=\mathrm{fl}(R_{ij}{x}_j)$ and then $s:=\mathrm{fl}(s+p)$; next compute $t:=\mathrm{fl}(c_i-s)$ and $x_i:=\mathrm{fl}(t\mathrm{/}{R}_{ii})$. If any product, addition, subtraction, or division fails, return the first range-failure. This is the descending recurrence of Forward and backward substitution for triangular linear systems with its rounded operation order made explicit.
   3. Compute ${\widehat{v}}_i$ for $i=0,\dots ,d-1$ as the rounded left-to-right dot product of row $i$ of the original $M$ with $x$: initialize the accumulator with $\mathrm{fl}(M_{i0}{x}_0)$ and for $j=1,\dots ,d-1$ compute $p:=\mathrm{fl}(M_{ij}{x}_j)$ and then replace the accumulator $s$ by $\mathrm{fl}(s+p)$. Set the resulting accumulator to ${\widehat{v}}_i$. Then compute each ${\widehat{r}}_i:=\mathrm{fl}({\widehat{v}}_i-b_i)$. Return the first range-failure if any operation fails.
   4. Request evaluator bounds $\bar{m}\ge \mathrm{\parallel }M{\mathrm{\parallel }}_F$, $\bar{x}\ge \mathrm{\parallel }x{\mathrm{\parallel }}_2$, $\bar{b}\ge \mathrm{\parallel }b{\mathrm{\parallel }}_2$, $\bar{v}\ge \mathrm{\parallel }\widehat{v}{\mathrm{\parallel }}_2$, and $\bar{r}\ge \mathrm{\parallel }\widehat{r}{\mathrm{\parallel }}_2$. If any request is unavailable, return (residual-unverified, x, null). Request a certified upper bound $\bar{W}$ for ${\gamma }_d\bar{m}\bar{x}+u(\bar{v}+\bar{b})$, where ${\gamma }_d:=du\mathrm{/}(1-du)$. If this request is unavailable, return (residual-unverified, x, null). Request a certified upper bound $\rho$ for $\bar{r}+\bar{W}$; if unavailable, return (residual-unverified, x, null). By A computed shifted-system residual has a certified upper bound, $\rho \ge \mathrm{\parallel }Mx-b{\mathrm{\parallel }}_2$. Return (success, x, rho) if $\rho \le ϵ$; otherwise return (residual-not-certified, x, rho).

This is a factor-only operation followed by an independently callable solve operation: the factor object accepts each later stored right-hand side and tolerance without receiving a precomputed list. Each solve performs one triangular solve and one residual check using the same stored factorization. A returned residual bound is an equation-residual certificate, not a forward-error estimate; converting it to forward error requires a separate inverse bound. No rank-revealing guarantee, sparse cost bound, or success guarantee for an ill-conditioned matrix is part of this procedure.

Factorization terminates with singular-factor, range-failure, factorization-failure, or a reusable factor object. Each call to the object's solve operation terminates with success, range-failure, residual-unverified, or residual-not-certified, together with any completed approximate solution and residual bound.

Use the one-stage map of The Rosenbrock–Euler step map with an exact or approximate Jacobian with the same $t,y,h,\gamma ,f$ but matrices $J$ and $\tilde{J}$. If $I_d-\gamma hJ$ and $I_d-\gamma h\tilde{J}$ are invertible, with inverse norms bounded by $K$ and $\tilde{K}$, then $\mathrm{\parallel }{\mathrm{\Phi }}_{\gamma ,h}^{\tilde{J}}(t,y)-{\mathrm{\Phi }}_{\gamma ,h}^J(t,y){\mathrm{\parallel }}_2\le \gamma h^{2}K\tilde{K}\mathrm{\parallel }\tilde{J}-J{\mathrm{\parallel }}_2\mathrm{\parallel }f(t,y){\mathrm{\parallel }}_2\mathrm{.}$ This is a one-step perturbation bound; it alone gives no uniform stability or convergence claim for a sequence of approximate Jacobians.

Let $d\ge 1$ and let $f:U\to {\mathbb{R}}^d$ be autonomous and $C^1$ on an open $U\subseteq {\mathbb{R}}^d$, let $y\in U$, and let $h>0$. Suppose $I_d-hDf(y)$ is invertible. Then the exact-arithmetic Rosenbrock–Euler step with $\gamma =1$ is the first Newton correction of the backward-Euler residual at the initial guess $y$. This does not assert that the corrected point solves the full backward-Euler equation.

Input: a dimension $d\ge 1$, an ODE field $f:D\to {\mathbb{R}}^d$ on the open domain $D\subseteq \mathbb{R}\times {\mathbb{R}}^d$, $(t_n,y_n)\in D$, $h>0$, $\gamma >0$, a finite real Jacobian matrix ${\tilde{J}}_n\in {\mathbb{R}}^{d\times d}$, a residual tolerance $\epsilon \ge 0$, a unit roundoff $u>0$ with $du<1$, the selected stored finite set $\mathbb{F}$ of Operationwise real and complex floating point arithmetic with rounded square roots and fused multiply-add, a terminating field-evaluation routine that returns either the exact value $f(t,y)$ or evaluation-failure for each requested $(t,y)\in D$, a terminating partial exact-storage check store, and a terminating certificate evaluator evaluator meeting the input contract of Reusable Householder QR solve for a shifted matrix with residual status. For an exact real scalar, vector, or matrix $v$, store(v) returns $v$ if every entry belongs to $\mathbb{F}$ and returns unavailable otherwise; it does not round. Apply it before each factorization call.

1. Form $M=I_d-\gamma h{\tilde{J}}_n$ by exact outer arithmetic and set $M_s:=\mathrm{store}(M)$. If storage returns unavailable, return $(y_n,\mathtt{\text{representation-unavailable}})$. Call factor(M_s,u,evaluator) of Reusable Householder QR solve for a shifted matrix with residual status. If it returns a factorization failure status, return $(y_n,\text{that status})$; otherwise retain the returned factor object.
2. Request the field evaluation at $(t_n,y_n)$. If it returns evaluation-failure, return $(y_n,\mathtt{\text{field-evaluation-failure}})$. Otherwise set $b=f(t_n,y_n)$ to the returned exact value and set $b_s:=\mathrm{store}(b)$. If storage returns unavailable, return $(y_n,\mathtt{\text{representation-unavailable}})$. Call the retained factor object's solve(b_s,\varepsilon). If the solve status is not success, return $(y_n,\text{that status})$ and retain its residual bound for diagnosis.
3. If the solve returns $k_n$ and a certified residual upper bound ${\rho }_n\le \epsilon$, return $(y_n+h{k}_n,\mathtt{\text{success}},{\rho }_n)$.

The name Rosenbrock–Euler applies when $\gamma =1$. The exact-arithmetic stage equation is $M{k}_n=f(t_n,y_n)$, exactly the map of The Rosenbrock–Euler step map with an exact or approximate Jacobian. Its residual $M{k}_n-f(t_n,y_n)$ is the update residual of the linear stage equation in the sense of Implicit residuals and their update and derivative units. For this mixed-arithmetic specification, evaluations and outer affine combinations are exact real operations and every matrix/RHS passed to factorization or solve must pass store; only the factorization, triangular solve, and residual dot products use rounded operations, while residual-bound formation is delegated to the terminating certificate evaluator. A full hardware-rounding theorem for outer operations needs additional certified evaluation and update-error terms. This implementation has one shifted-matrix factorization and one stage solve. It does not turn a small residual into a forward-error claim without an inverse bound.

Input: a dimension $d\ge 1$, an ODE field $f:D\to {\mathbb{R}}^d$ on the open domain $D\subseteq \mathbb{R}\times {\mathbb{R}}^d$, $(t_n,y_n)\in D$, $h>0$, a finite real state-Jacobian approximation ${\tilde{J}}_n\in {\mathbb{R}}^{d\times d}$, a finite approximation ${\tilde{T}}_n\in {\mathbb{R}}^d$ to ${\mathrm{\partial }}_{t}f(t_n,y_n)$, stage residual tolerances ${\epsilon }_1,{\epsilon }_2\ge 0$, a unit roundoff $u>0$ with $du<1$, the selected stored finite set $\mathbb{F}$ of Operationwise real and complex floating point arithmetic with rounded square roots and fused multiply-add, a terminating field-evaluation routine evaluate_f that returns either the exact value $f(t,y)$ or evaluation-failure for each requested $(t,y)\in D$, a terminating domain-membership routine in_domain returning true exactly for points in $D$ and false otherwise, a terminating store routine, and terminating certificate evaluator evaluator meeting the input contract of Reusable Householder QR solve for a shifted matrix with residual status. The exact-storage check returns $v$ if every entry is in $\mathbb{F}$ and unavailable otherwise; it does not round. Put $\alpha :=\frac{1}{2+\sqrt{2}}=1-\frac{1}{\sqrt{2}},W:=I_d-\alpha h{\tilde{J}}_n\mathrm{.}$ The two-stage method is $W{k}_1=f(t_n,y_n)+\alpha h{\tilde{T}}_n,$ $F_1=f\text{\hspace{0.17em}⁣}(t_n+\frac{h}{2},y_n+\frac{h}{2}{k}_1),W{q}_2=F_1-k_1,k_2=q_2+k_1,$ $y_{n+1}=y_n+h{k}_2\mathrm{.}$

The residuals $W{k}_1-[f(t_n,y_n)+\alpha h{\tilde{T}}_n]$ and $W{q}_2-(F_1-k_1)$ are update residuals of their linear stage equations in the sense of Implicit residuals and their update and derivative units. The coefficient $\alpha$ and the autonomous two-stage formula are the second-order Rosenbrock formula recorded for ode23s in WIAS Example 2.37. The first-stage $\alpha h{\tilde{T}}_n$ term specifies the nonautonomous extension used here; its order is independently established under the stated boundedness hypotheses.

1. Form $W$ by exact outer arithmetic and set $W_s:=\mathrm{store}(W)$. If storage returns unavailable, return $(y_n,\mathtt{\text{representation-unavailable}})$. Call factor(W_s,u,evaluator) of Reusable Householder QR solve for a shifted matrix with residual status. If it returns a factorization failure status, return $(y_n,\mathtt{\text{factorization-failure}})$; otherwise retain its factor object.
2. Call evaluate_f(t_n,y_n). If it returns evaluation-failure, return $(y_n,\mathtt{\text{input-or-representation-failure}})$. Otherwise denote its returned exact value by $f_0$, form $b_1=f_0+\alpha h{\tilde{T}}_n$ by exact outer arithmetic, and set $b_{1,s}:=\mathrm{store}(b_1)$. If storage returns unavailable, return $(y_n,\mathtt{\text{input-or-representation-failure}})$. Call the retained factor object's solve(b_{1,s},\varepsilon_1). If it is not certified, return $(y_n,\text{stage-1 status})$. Otherwise set $k_1$ to its returned vector and ${\rho }_1$ to its certified residual bound.
3. Call in_domain on $(t_n+h\mathrm{/}2,y_n+h{k}_1\mathrm{/}2)$. If it returns false, return $(y_n,\mathtt{\text{invalid-stage-domain}})$. Otherwise call evaluate_f(t_n+h/2,y_n+hk_1/2). If it returns evaluation-failure, return $(y_n,\mathtt{\text{stage-2-input-or-representation-failure}})$. Otherwise set $F_1$ to the returned exact value, form $b_2=F_1-k_1$ exactly, and set $b_{2,s}:=\mathrm{store}(b_2)$. If storage returns unavailable, return $(y_n,\mathtt{\text{stage-2-input-or-representation-failure}})$. Call the same factor object's solve(b_{2,s},\varepsilon_2). If this solve is not certified, return $(y_n,\text{stage-2 status})$. Otherwise set $q_2$ to its returned vector and ${\rho }_2$ to its certified residual bound.
4. Set $k_2=q_2+k_1$ and return $(y_n+h{k}_2,\mathtt{\text{success}},{\rho }_1,{\rho }_2)$ with the two certified update-residual bounds.

Where all the indicated state and time partial derivatives exist, the exact-J branch takes ${\tilde{J}}_n=D_{y}f(t_n,y_n)$ and ${\tilde{T}}_n={\mathrm{\partial }}_{t}f(t_n,y_n)$. The W-method branch permits any bounded matrix approximations satisfying the hypotheses of the local and global order results; it receives no scalar A/L-stability guarantee from the exact-J result. The factor matrix is the same in both stages, so the procedure performs one shifted-matrix factorization and two stage solves. For this mixed-arithmetic specification, field evaluations and outer affine combinations are exact real operations, and every matrix/RHS passed to factorization or solve must pass store; factorization, triangular solves, and residual dot products use the operationwise model, while residual-bound formation is delegated to the terminating certificate evaluator. A full hardware-rounding theorem for outer operations would need additional certified evaluation and update-error terms. A finite arithmetic residual is not a forward-error bound without an inverse estimate.

Let $D\subseteq \mathbb{R}\times {\mathbb{R}}^d$ be open and $f:D\to {\mathbb{R}}^d$ be $C^2$ (componentwise, as in $C^k$ Euclidean maps and diffeomorphisms). Let $y(s)$ solve $y^{\mathrm{\prime }}=f(s,y)$ on $[t,t+h]$, $h>0$, and set $J=D_{y}f(t,y(t))$. Suppose the trajectory lies in a region where $\mathrm{\parallel }f{\mathrm{\parallel }}_2\le F$, $\mathrm{\parallel }{D}_{y}f{\mathrm{\parallel }}_2\le L$, and $\mathrm{\parallel }{f}_t+D_{y}ff{\mathrm{\parallel }}_2\le M_1$. Fix $\gamma >0$ and a real matrix $\tilde{J}$. Suppose $A=I-\gamma hJ$ and $\tilde{A}=I-\gamma h\tilde{J}$ are invertible with $\mathrm{\parallel }{A}^{-1}{\mathrm{\parallel }}_2\le {\kappa }_0$ and $\mathrm{\parallel }{\tilde{A}}^{-1}{\mathrm{\parallel }}_2\le {\kappa }_1$. Write $\delta =\mathrm{\parallel }\tilde{J}-J{\mathrm{\parallel }}_2$.

For the exact-arithmetic Rosenbrock stage map ${\mathrm{\Phi }}_{\gamma ,h}^{\tilde{J}}$, and a computed stage $\widehat{k}$ whose update residual obeys $\mathrm{\parallel }\tilde{A}\widehat{k}-f(t,y(t)){\mathrm{\parallel }}_2\le \rho$, the update ${\widehat{y}}_{+}=y(t)+h\widehat{k}$ satisfies $\mathrm{\parallel }{\widehat{y}}_{+}-y(t+h){\mathrm{\parallel }}_2\le C_E{h}^2+{\kappa }_{1}h\rho ,C_E:=\gamma {\kappa }_{0}LF+\frac{M_1}{2}+\gamma {\kappa }_0{\kappa }_1\delta F\mathrm{.}$ The last term is the separate Jacobian-approximation allowance. The residual is an equation residual in the units of Implicit residuals and their update and derivative units, not a forward-error estimate by itself.

Let $D\subseteq \mathbb{R}\times {\mathbb{R}}^d$ be open, $f:D\to {\mathbb{R}}^d$ be $C^2$, and $y(s)$ solve $y^{\mathrm{\prime }}=f(s,y)$ on $[t,t+h]$. Assume a compact convex tube $Q\subset D$ contains the trajectory, the stage midpoint segments, and all vectors used below. On $Q$, suppose

where $D^{2}f$ is the second total derivative on the joint variable $(t,y)$, and $\mathrm{\parallel }{D}^{2}f{\mathrm{\parallel }}_{\mathrm{o}\mathrm{p}}:={\sup }_{\mathrm{\parallel }u\mathrm{\parallel }=\mathrm{\parallel }v\mathrm{\parallel }=1}\mathrm{\parallel }{D}^{2}f[u,v]{\mathrm{\parallel }}_2$ for the Euclidean product norm. Fix $\alpha =1-1\mathrm{/}\sqrt{2}$, and suppose $\mathrm{\parallel }\tilde{J}{\mathrm{\parallel }}_2\le B_J$, $\mathrm{\parallel }\tilde{T}{\mathrm{\parallel }}_2\le B_T$, and $\mathrm{\parallel }\alpha h\tilde{J}{\mathrm{\parallel }}_2\le q<1$. Put $\kappa =(1-q{)}^{-1}$ and $M=T_0+LF$. The shifted inverse has norm at most $\kappa$ by A small shifted Jacobian has an explicit inverse bound.

For exact-arithmetic stages $k_1,q_2,k_2$ of Two-stage second-order Rosenbrock W-method with shared factorization (with its exact-storage checks successful), define $C_1=B_{J}F+B_T,K_1=F+\alpha h\kappa C_1,C_F=\frac{\alpha \kappa L{C}_1}{2}+\frac{H(1+K_1^2)}{8}\mathrm{.}$ $C_K={\alpha }^2{\kappa }^2{B}_J{C}_1+\frac{\alpha \kappa B_{J}M}{2}+\kappa C_F,M_2=H(1+F^2)+L{T}_0+L^{2}F\mathrm{.}$ Then the exact one-step defect is at most

If the computed first and second shifted solves have certified update-residual bounds ${\rho }_1$ and ${\rho }_2$, then the computed update satisfies $\mathrm{\parallel }{\widehat{y}}_{+}-y(t+h){\mathrm{\parallel }}_2\le (C_K+M_2\mathrm{/}6)h^3+C_r{h}^2{\rho }_1+\kappa h{\rho }_2,C_r:={\kappa }^2(L\mathrm{/}2+\alpha B_J)\mathrm{.}$ Thus ${\rho }_1=O(h)$ and ${\rho }_2=O(h^2)$ preserve a cubic local defect. These are equation-residual conditions, not stage forward-error claims without the inverse bound.

At any base point in the tube, the same two-stage residual perturbation argument gives the solver-to-exact-map bound $\mathrm{\parallel }{\widehat{\mathrm{\Phi }}}_h(t,y)-{\mathrm{\Phi }}_h(t,y){\mathrm{\parallel }}_2\le C_r{h}^2{\rho }_1+\kappa h{\rho }_2$, whenever both midpoint evaluations lie in $Q$.