Implicit stiff time stepping and algebraic error — Examples

For all $x,y\in \mathbb{R}$, $(f(x)-f(y))(x-y)=-(x-y{)}^2(x^2+xy+y^2)\le 0,$ because $x^2+xy+y^2=(x+y\mathrm{/}2{)}^2+3y^2\mathrm{/}4\ge 0$. This polynomial field is continuous on $\mathbb{R}$ and is globally dissipative with $\mu =0$.

By [F1], the equation $z=1-z^3$, equivalently $z+z^3=1$, has exactly one real root. At the supplied initial iterate $x_0=0$, [F3] gives $G(0)=0-1-(-0)=-1$, so $\mathrm{\parallel }G(0)\mathrm{\parallel }=1>0.1$.

Since the current residual exceeds tolerance and the correction budget is already zero, the finite solver returns iteration-budget-exhausted, despite the unique BE root proved in step 2.1. No forward-error statement follows from this status.

Here the test argument is $z=-\lambda h=-100$. The trapezoidal rule equation $y_{n+1}=y_n+\frac{h}{2}(-\lambda y_n-\lambda y_{n+1})$ rearranges to the factor $R_{\mathrm{t}\mathrm{r}\mathrm{a}\mathrm{p}}(-100)=\frac{1-100\mathrm{/}2}{1+100\mathrm{/}2}=-\frac{49}{51}\mathrm{.}$ Thus, for $y_0\mathrm{\ne }0$, successive trapezoidal values alternate sign, while the magnitude after each step is multiplied by $49\mathrm{/}51\approx 0.961$; for $y_0=0$, every iterate remains zero.

Backward Euler gives $R_{\mathrm{B}\mathrm{E}}(-100)=1\mathrm{/}101\approx 0.00990$. The proved scalar functions give $R_S(-100)=\frac{1-100(1-2\gamma )}{(1+100\gamma {)}^2}\approx -0.0441,R_R(-100)=\frac{1-100\mathrm{/}3}{1+200\mathrm{/}3+10000\mathrm{/}6}\approx -\mathrm{0.0186.}$ Each L-stable method damps this stiff mode strongly at this step size; the trapezoidal factor instead changes sign and has magnitude near one. These are scalar linear-test calculations only and imply no nonlinear stability comparison.

The scalar field is strongly dissipative with $\mu =1$, since for all $x,y$, $(f(x)-f(y))(x-y)=-(x-y{)}^2(1+x^2+xy+y^2)\le -(x-y{)}^2,$ where $x^2+xy+y^2\ge 0$. Hence [F1] gives a unique BE root on all of $\mathbb{R}$. The root equation is $z+2z+2z^3=1$, or $3z+2z^3=1$.

At ${\tilde{u}}_1=1\mathrm{/}3$, the left side of the root equation is $1+2\mathrm{/}27$, so the update residual is ${\eta }_0=2\mathrm{/}27$ and the derivative residual is ${\eta }_0\mathrm{/}h=1\mathrm{/}27$. Here $q=1\mathrm{/}(1+2)=1\mathrm{/}3$. Applying [F2] gives $\mathrm{\mid }{\tilde{u}}_1-u_1\mathrm{\mid }\le q\mathrm{\mid }{\eta }_0\mathrm{\mid }=\frac{1}{3}\frac{2}{27}=\frac{2}{81}\mathrm{.}$

Differentiating the polynomial components gives the displayed $Df$ by [F5]. For an increment $v$, the remainder $f(y+v)-f(y)-Df(y)v$ is $({10}^4{v}_2{v}_3,-{10}^4{v}_2{v}_3-3\cdot {10}^7{v}_2^2,3\cdot {10}^7{v}_2^2{)}^T$, whose Euclidean norm is bounded by $C\mathrm{\parallel }v{\mathrm{\parallel }}_2^2$; hence $Df(y)$ is the total derivative and its affine entries are continuous. Thus $J_G=I-hDf$; by [F4], both $f$ and $G(y)=y-y_0-hf(y)$ are $C^1$, as required by [F2]. Summing the components of $f$ cancels all terms, so a root satisfies ${\sum }_i(y_i-y_{0,i})=0$ and preserves total concentration.

At the start, $f(y_0)=(-0.04,0.04,0{)}^T$, so $G(y_0)=(0.00004,-0.00004,0{)}^T$ and $\mathrm{\parallel }G(y_0){\mathrm{\parallel }}_2=\sqrt{2}(4\cdot {10}^{-5})=5.65685425\cdot {10}^{-5}>{10}^{-12}$ by [F6]. Thus neither mode returns at count zero.

Applying the stated Newton correction formula with the displayed Jacobian and binary64 partial-pivoting elimination gives the Newton residual sequence in the Numerical record. It first falls below ${10}^{-12}$ at count five; the residuals at counts zero through four exceed the tolerance. Since the budget is 12, this binary64 implementation returns root after five corrections. Each correction refreshes and factors one Jacobian, then solves once, giving five assemblies/factorizations and five solves.

With chord refresh period eight, [F2] refreshes at count zero and retains that factorization for corrections zero through seven. The residuals at counts zero through eight all exceed ${10}^{-12}$, and the budget is eight, so this binary64 implementation returns iteration-budget-exhausted. This uses one assembly/factorization and eight linear solves. The mode is unconstrained, so negative intermediate second components do not violate its domain ${\mathbb{R}}^3$.