Two-stage second-order Rosenbrock W-method with shared factorization

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.