Implicit stiff time stepping and algebraic error

Let $G:U\to {\mathbb{R}}^m$ be the residual map of an implicit equation, where $U\subseteq {\mathbb{R}}^m$ is the admissible set for the unknown vector. An exact solution is a point $z^{\ast }\in U$ with $G(z^{\ast })=0$.

An equation residual need not have the same units as its unknown. When the equation is normalized so that each component of $G(z)$ has the units of the corresponding component of $z$, the update residual at an approximate solution $\tilde{z}\in U$ is $G(\tilde{z})$. For example, $z^2-a$ has squared state units when $a$ does, so it must first be normalized to use this update-residual convention. If the original equation was divided by a positive step scale $\sigma$, then the derivative residual is $G(\tilde{z})\mathrm{/}\sigma$, and the corresponding update residual is $\sigma$ times that derivative residual. A residual norm is not a forward-error bound unless a separate inverse or stability estimate converts it to $\mathrm{\parallel }\tilde{z}-z^{\ast }\mathrm{\parallel }$.

For a first-order vector field $f:D\to {\mathbb{R}}^d$ on an open time-state domain $D\subseteq \mathbb{R}\times {\mathbb{R}}^d$ (First-order systems, initial value problems, and solutions on intervals), take a time $t$, input $y\in {\mathbb{R}}^d$, and step $h>0$ with $(t,y)\in D$. For backward Euler, an endpoint $z$ is admissible only when $(t+h,z)\in D$, and its update residual is

$G_{\mathrm{B}\mathrm{E}}(z)=z-y-hf(t+h,z),$

so the derivative residual is $G_{\mathrm{B}\mathrm{E}}(z)\mathrm{/}h$.

For an $s$-stage implicit Runge–Kutta tableau $(A,b,c)$, a stage vector $Y=(Y_1,\dots ,Y_s)$ is admissible only when $(t+c_{j}h,Y_j)\in D$ for every $j$. On that admissible set, its stage update residual has blocks

$G_i(Y)=Y_i-y-h{\sum }_{j=1}^s{a}_{ij}f(t+c_{j}h,Y_j),i=1,\dots ,s\mathrm{.}$

The stage derivative residual is $G(Y)\mathrm{/}h$. These conventions distinguish an equation residual from the distance to its exact root.

Let $d\ge 1$, $\kappa >0$, and let $H:{\mathbb{R}}^d\to {\mathbb{R}}^d$ be continuous. Suppose

$⟨H(x)-H(y),x-y⟩\ge \kappa \mathrm{\parallel }x-y{\mathrm{\parallel }}_2^2(x,y\in {\mathbb{R}}^d)\mathrm{.}$

Then $H$ is bijective. Its inverse satisfies

$\mathrm{\parallel }{H}^{-1}(u)-H^{-1}(v){\mathrm{\parallel }}_2\le {\kappa }^{-1}\mathrm{\parallel }u-v{\mathrm{\parallel }}_2(u,v\in {\mathbb{R}}^d)\mathrm{.}$

Let $f:D\to {\mathbb{R}}^d$ define a first-order initial-value problem on an open time-state domain $D\subseteq \mathbb{R}\times {\mathbb{R}}^d$. Given a time $t_n$, step $h_n>0$, and state $y_n$ with $(t_n,y_n)\in D$, an admissible endpoint $y_{n+1}$ must satisfy $(t_n+h_n,y_{n+1})\in D$. The backward Euler endpoint is any such $y_{n+1}$ satisfying

$y_{n+1}=y_n+h_{n}f(t_n+h_n,y_{n+1})\mathrm{.}$

Its update residual at an approximate endpoint $z$ with $(t_n+h_n,z)\in D$ is

$G_n^{\mathrm{B}\mathrm{E}}(z)=z-y_n-h_{n}f(t_n+h_n,z)\mathrm{.}$

When $f$ is defined at the endpoint, the derivative residual is $G_n^{\mathrm{B}\mathrm{E}}(z)\mathrm{/}{h}_n$. The equation can have no, one, or several roots on a restricted admissible branch set; this definition by itself asserts neither existence nor uniqueness.

Let $d,s\ge 1$, $A=(a_{ij})\in {\mathbb{R}}^{s\times s}$, $b,c\in {\mathbb{R}}^s$, and let $f:D\to {\mathbb{R}}^d$ be a vector field on an open time-state domain $D\subseteq \mathbb{R}\times {\mathbb{R}}^d$ (First-order systems, initial value problems, and solutions on intervals). For an initial state $y\in {\mathbb{R}}^d$ and time $t$ with $(t,y)\in D$, and a step $h>0$, a stage vector $Y=(Y_1,\dots ,Y_s)\in ({\mathbb{R}}^d{)}^s$ is admissible only when $(t+c_{j}h,Y_j)\in D$ for every $j$. An admissible vector solves

$Y_i=y+h{\sum }_{j=1}^s{a}_{ij}f(t+c_{j}h,Y_j),i=1,\dots ,s\mathrm{.}$

When a solution $Y$ is selected, the method update is

${\mathrm{\Psi }}_h(t,y)=y+h{\sum }_{i=1}^s{b}_{i}f(t+c_{i}h,Y_i)\mathrm{.}$

On the admissible stage set, the stage update residual is the vector with blocks

$G_i(Y)=Y_i-y-h{\sum }_{j=1}^s{a}_{ij}f(t+c_{j}h,Y_j)\mathrm{.}$

The associated stage derivative residual is $G(Y)\mathrm{/}h$. The tableau is stiffly accurate when $b_j=a_{sj}$ for all $j$ and $c_s=1$; in that case the last stage $Y_s$ equals the method update. The equations may have no, one, or several solutions in a proposed branch set; this definition alone asserts no solvability.

Let $\gamma =1-1\mathrm{/}\sqrt{2}$. The selected two-stage singly diagonally implicit Runge–Kutta method has

$A=\left(\begin{array}{cc}\gamma & 0\\ 1-\gamma & \gamma \end{array}\right),b=\left(\begin{array}{c}1-\gamma \\ \gamma \end{array}\right),c=\left(\begin{array}{c}\gamma \\ 1\end{array}\right)\mathrm{.}$

Applied to $y^{\mathrm{\prime }}=f(t,y)$, its stages and update are those of Implicit Runge–Kutta stage equations and update. Both diagonal entries of $A$ equal $\gamma$, the row sums are $c$, and the final row of $A$ is $b^{\mathsf{T}}$.

The two-stage Radau IIA method has

$A=\left(\begin{array}{cc}5\mathrm{/}12& -1\mathrm{/}12\\ 3\mathrm{/}4& 1\mathrm{/}4\end{array}\right),b=\left(\begin{array}{c}3\mathrm{/}4\\ 1\mathrm{/}4\end{array}\right),c=\left(\begin{array}{c}1\mathrm{/}3\\ 1\end{array}\right)\mathrm{.}$

Applied to $y^{\mathrm{\prime }}=f(t,y)$, its stages and update are those of Implicit Runge–Kutta stage equations and update. The row sums of $A$ equal $c$, and the final row of $A$ is $b^{\mathsf{T}}$.

For a Runge–Kutta tableau $(A,b,c)$, apply its stage equations to the scalar test equation $y^{\mathrm{\prime }}=\lambda y$, with $z=h\lambda$. Whenever $I-zA$ is invertible, the stages satisfy $Y=(I-zA{)}^{-1}\mathbf{1}y$, and the scalar stability function is

$R(z)=1+z{b}^{\mathsf{T}}(I-zA{)}^{-1}\mathbf{1.}$

The scalar stability region is the set of $z\in \mathbb{C}$ for which the scalar update has $\mathrm{\mid }R(z)\mathrm{\mid }\le 1$. The method is A-stable when its rational stability function has no pole in $\{\mathrm{Re}z\le 0\}$ and $\mathrm{\mid }R(z)\mathrm{\mid }\le 1$ throughout that closed half-plane. It is L-stable when it is A-stable and $R(z)\to 0$ as $\mathrm{\mid }z\mathrm{\mid }\to \mathrm{\infty }$ with $\mathrm{Re}z\le 0$.

These are properties of the scalar linear test equation and the tableau's rational function. They do not assert nonlinear contractivity or control of nonnormal matrix growth.

Let $m,d\ge 1$, $h>0$, $A\in {\mathbb{R}}^{m\times m}$, and centers $c=(c_1,\dots ,c_m)\in ({\mathbb{R}}^d{)}^m$. On block vectors define

$\mathrm{\parallel }Y{\mathrm{\parallel }}_{\mathrm{\infty },2}:={\max }_{1\le i\le m}\mathrm{\parallel }{Y}_i{\mathrm{\parallel }}_2,\mathrm{\parallel }A{\mathrm{\parallel }}_{\mathrm{\infty }}:={\max }_i{\sum }_{j=1}^m\mathrm{\mid }{a}_{ij}\mathrm{\mid }\mathrm{.}$

Let $F:\mathcal{B}\to ({\mathbb{R}}^d{)}^m$, with components $F_i:\mathcal{B}\to {\mathbb{R}}^d$, be defined on the closed product ball $\mathcal{B}=\{Y:\mathrm{\parallel }Y-c{\mathrm{\parallel }}_{\mathrm{\infty },2}\le R\}$, where $R>0$. Suppose constants $M,L\ge 0$ satisfy, for all $Y,Z\in \mathcal{B}$,

$\mathrm{\parallel }F(Y){\mathrm{\parallel }}_{\mathrm{\infty },2}\le M,\mathrm{\parallel }F(Y)-F(Z){\mathrm{\parallel }}_{\mathrm{\infty },2}\le L\mathrm{\parallel }Y-Z{\mathrm{\parallel }}_{\mathrm{\infty },2},$

and $h\mathrm{\parallel }A{\mathrm{\parallel }}_{\mathrm{\infty }}M\le R$, $q:=h\mathrm{\parallel }A{\mathrm{\parallel }}_{\mathrm{\infty }}L<1$. Define $\mathrm{\Phi }(Y)=c+h(A\otimes I_d)F(Y)$ and $G(Y)=Y-\mathrm{\Phi }(Y)$. Then $\mathrm{\Phi }$ has a unique fixed point $Y^{\ast }\in \mathcal{B}$, equivalently $G(Y^{\ast })=0$. Every approximate stage vector $\tilde{Y}\in \mathcal{B}$ satisfies

$\mathrm{\parallel }\tilde{Y}-Y^{\ast }{\mathrm{\parallel }}_{\mathrm{\infty },2}\le \frac{\mathrm{\parallel }G(\tilde{Y}){\mathrm{\parallel }}_{\mathrm{\infty },2}}{1-q}\mathrm{.}$

If the derivative-scaled stage residual is $r=G(\tilde{Y})\mathrm{/}h$, the right side is $h\mathrm{\parallel }r{\mathrm{\parallel }}_{\mathrm{\infty },2}\mathrm{/}(1-q)$. The hypotheses $h\mathrm{\parallel }A{\mathrm{\parallel }}_{\mathrm{\infty }}L<1$ and the invariant-ball bound are local solvability conditions; this result gives no global branch outside $\mathcal{B}$.

Let $m\ge 1$, let $U\subseteq {\mathbb{R}}^m$ be open, and let $G:U\to {\mathbb{R}}^m$ be $C^1$, $z_0\in U$, residual tolerance $\epsilon \ge 0$, correction budget $N\in {\mathbb{N}}_0$, and refresh period $q\in \mathbb{N}$, $q\ge 1$. Use the coordinate Jacobian convention of The Jacobian matrix of partial derivatives and the gradient in the scalar-valued case. Since $G$ is $C^1$, A total derivative computes every directional derivative, and its matrix is the Jacobian identifies this Jacobian with the matrix of $DG(z)$; write it as $J(z)$. Residual and Jacobian evaluations, norm comparisons and domain tests are returning primitives; no bound on their internal cost is specified. In chord mode also supply a terminating exact linear-system backend: on each square Jacobian it reports singularity if the matrix is singular, and otherwise returns a retained factor object, and each solve on that object terminates with the exact unique solution for its right-hand side. Factorization and solve work are charged individually. The output is a last valid iterate $z\in U$, a status, and the actual residual norm $\mathrm{\parallel }G(z){\mathrm{\parallel }}_2$.

In updated-Newton mode, invoke Finite-budget Newton for a square vector residual with $G,z_0,\epsilon ,N$. It recomputes the Jacobian at every correction and returns its stated residual/budget/failure status. Re-evaluate $G$ at the returned iterate to provide the output residual norm. Count each Jacobian assembly, factorization attempt, and linear solve actually performed.

In chord mode, for $k=0,1,\dots ,N$, first evaluate $G(z_k)$ and put $r_k=\mathrm{\parallel }G(z_k){\mathrm{\parallel }}_2$. If $r_k\le \epsilon$, return $(z_k,\mathtt{\text{root}},r_k)$; this status means only that the residual tolerance was met. If $k=N$, return $(z_k,\mathtt{\text{iteration-budget-exhausted}},r_k)$. Otherwise, at $k=0,q,2q,\dots$, refresh by forming $J_k=DG(z_k)$, attempting to factor it, and retaining the factorization until the next scheduled refresh. The factorization backend must either report singularity or return a factor object whose solve operation gives the exact unique solution of $J_{k}s=-G(z_k)$. If it reports singularity, return $(z_k,\mathtt{\text{singular-jacobian}},r_k)$. At each correction obtain $s_k$ by calling the solve operation on that retained factor object with right-hand side $-G(z_k)$; do not refactor or invoke a fresh matrix-only solve. The exact-solve backend contract makes this the mathematical correction in Newton correction for a square vector residual. If the correction $s_k$ has $z_k+s_k\mathrm{\notin }U$, return $(z_k,\mathtt{\text{invalid-trial}},r_k)$; otherwise set $z_{k+1}=z_k+s_k$ and continue. A refresh period $q=1$ updates the Jacobian every correction; $q>1$ is a chord iteration.

Both modes make no claim that a finite budget reaches an exact root or that a small residual is a forward-error bound. Record the number of Jacobian assemblies, factorization attempts, and linear solves actually performed; a dense LU correction uses one forward and one backward triangular solve. Sparse cost depends on the actual stage matrix and factorization backend and is not bounded here.

Let $d\ge 1$, let $f:\mathbb{R}\times {\mathbb{R}}^d\to {\mathbb{R}}^d$ be continuous, and suppose for every fixed $t$ and all $x,y\in {\mathbb{R}}^d$,

$⟨f(t,x)-f(t,y),x-y⟩\le -\mu \mathrm{\parallel }x-y{\mathrm{\parallel }}_2^2$

for some $\mu \ge 0$. For every $h>0$, define $H_{t,h}(x)=x-hf(t,x)$. Then $H_{t,h}$ is bijective and its inverse $R_{t,h}=H_{t,h}^{-1}$, the backward Euler step map, satisfies

$\mathrm{\parallel }{R}_{t,h}(u)-R_{t,h}(v){\mathrm{\parallel }}_2\le \frac{1}{1+h\mu }\mathrm{\parallel }u-v{\mathrm{\parallel }}_2\mathrm{.}$

In particular, for each input state $u$, the equation $z=u+hf(t,z)$ has exactly one solution on all of ${\mathbb{R}}^d$. No $hL<1$ contraction assumption is needed for this dissipative branch.

Let $\gamma =1-1\mathrm{/}\sqrt{2}$, and use the tableau in Selected two-stage singly diagonally implicit Runge–Kutta tableau. The coefficient satisfies $0<\gamma <1$. For a $C^2$ vector field on a neighborhood of the exact solution segment, the local small-step stage branch specified by Small-step implicit Runge–Kutta stages have a unique branch and controlled residual error has local defect $O(h^3)$; hence the method has classical order at least two. Its scalar stability function is $R(z)=\frac{1+(1-2\gamma )z}{(1-\gamma z{)}^2},$ and it is L-stable in the scalar sense of Scalar stability function, A-stability, and L-stability for Runge–Kutta methods. This is only a scalar linear-test-equation assertion, not a nonlinear contractivity claim.

For the tableau in Two-stage Radau IIA tableau, on a $C^3$ vector field and the local small-step branch specified by Small-step implicit Runge–Kutta stages have a unique branch and controlled residual error, the local defect is $O(h^4)$, so the method has classical order at least three. Its scalar stability function is $R(z)=\frac{1+z\mathrm{/}3}{1-2z\mathrm{/}3+z^2\mathrm{/}6},$ and it is L-stable in the scalar sense of Scalar stability function, A-stability, and L-stability for Runge–Kutta methods. No nonlinear contractivity assertion follows.

Assume the hypotheses of Backward Euler has a globally defined nonexpansive resolvent (contractive when strongly dissipative) with a common $\mu \ge 0$, let $N\in {\mathbb{N}}_0$, and use a constant $h>0$ grid $t_n=t_0+nh$ for $0\le n\le N$. Let $u_0\in {\mathbb{R}}^d$ and let $u_{n+1}$ be the exact backward Euler sequence, $u_{n+1}=R_{t_{n+1},h}(u_n)$. For any approximate sequence ${\tilde{u}}_n\in {\mathbb{R}}^d$, define its update residual ${\eta }_n={\tilde{u}}_{n+1}-{\tilde{u}}_n-hf(t_{n+1},{\tilde{u}}_{n+1})\mathrm{.}$ With $q=(1+h\mu {)}^{-1}\le 1$, for $0\le n\le N$, $\mathrm{\parallel }{\tilde{u}}_N-u_N{\mathrm{\parallel }}_2\le q^N\mathrm{\parallel }{\tilde{u}}_0-u_0{\mathrm{\parallel }}_2+{\sum }_{j=0}^{N-1}{q}^{N-j}\mathrm{\parallel }{\eta }_j{\mathrm{\parallel }}_2\mathrm{.}$ This bound uses the unscaled update residual. If $r_n={\eta }_n\mathrm{/}h$ is the derivative residual, its contribution is $h{q}^{N-j}\mathrm{\parallel }{r}_j{\mathrm{\parallel }}_2$. The result compares to the exact discrete BE sequence, not directly to the continuous solution.