Rosenbrock and linearly implicit methods — Examples

Define the altered variant as stated. Since $W=1-\alpha h$, its first stage is $k_1=y\mathrm{/}(1-\alpha h)$. The midpoint value is $F_1=y+(h\mathrm{/}2)k_1$, and the altered second solve gives ${\mathrm{\ell }}_2=F_1\mathrm{/}(1-\alpha h)$.

Define $Y(s)=\exp (s)$ for $s\in \mathbb{R}$. By [F4], $Y^{\mathrm{\prime }}=Y$, and [F6] gives $Y(0)=1$; thus this is an exact solution of the IVP (unique on $\mathbb{R}$ by [F3]).

Put $r=(1-\alpha h{)}^{-1}$. The exact identity $(1-x{)}^{-1}=1+x+x^2+x^3\mathrm{/}(1-x)$ shows $r=1+\alpha h+{\alpha }^2{h}^2+O(h^3)$; for $0<h\le 1\mathrm{/}(2\alpha )$ the remainder is bounded by $2{\alpha }^3{h}^3$. Squaring this finite expansion gives $r^2=1+2\alpha h+3{\alpha }^2{h}^2+O(h^3)$, with a bounded remainder on the same interval. Since $F_1=y+(h\mathrm{/}2)k_1=y(1+hr\mathrm{/}2)$, the altered second stage is ${\mathrm{\ell }}_2=F_{1}r=y(r+h{r}^2\mathrm{/}2)=y(1+(\alpha +1\mathrm{/}2)h+({\alpha }^2+\alpha )h^2+O(h^3))$. Multiplying by $h$ and adding $y$ gives the displayed one-step expansion.

By [F4], $Y^{\mathrm{\prime }\mathrm{\prime }\mathrm{\prime }}(s)=\exp (s)$; since $0\le s\le h\le 1$, [F5] gives $Y^{\mathrm{\prime }\mathrm{\prime }\mathrm{\prime }}(s)\le \exp (1)$.

The derivatives of $Y$ are continuous by [F4] and integrable on each compact subinterval by [F7]. Three applications of [F2] give $\exp (h)=1+h+h^2\mathrm{/}2+R$, where $R={\int }_0^h(h-s{)}^2{Y}^{\mathrm{\prime }\mathrm{\prime }\mathrm{\prime }}(s)\mathrm{/}2ds$ and $\mathrm{\mid }R\mathrm{\mid }\le \exp (1)h^3\mathrm{/}6$. Subtracting this expansion from the altered update in step 2.1 leaves $\alpha h^2+O(h^3)$, which is not bounded by $C{h}^3$ for any fixed $C$ as $h\downarrow 0$ because $\alpha >0$. Thus the altered scheme fails the cubic local-error condition for order two.

The proper method's second right-hand side and corrected output are exactly the two terms changed in this variant by [F1]. Step 3.1 proves the altered formula has no cubic local defect, and this item makes no global-rate claim for it.

Use the autonomous extension $f(t,y)=Jy$ on $D=\mathbb{R}\times {\mathbb{R}}^2$, with in_domain returning true and exact field evaluation and storage checks. The W approximation is $\tilde{J}=0$, so $W=I_2$. For the first active column $x=(1,0{)}^T$, [F3] has $s=1$, $\widehat{r}=1$, $\widehat{\alpha }=-1$, $v=(2,0{)}^T$, $w=4$, and $\widehat{\tau }=1\mathrm{/}2$. Its trailing-block dot is zero, so the output block is $\mathrm{diag}(-1,1)$. The second one-dimensional active column similarly has $v=2,w=4,\widehat{\tau }=1\mathrm{/}2$ and produces $R=-I_2$. By [F5] every exact representable intermediate is unchanged; all nonzero intermediates are among the listed stored normal values, no range branch fails, and both pivots are nonzero. Thus the one factorization succeeds.

The first RHS is $f(y)=(-1,-4{)}^T$. In the first reflector replay the dot is $-2$ and the active vector $p=\widehat{\tau }v$ is $(1,0{)}^T$, so the update sends $(-1,-4{)}^T$ to $(1,-4{)}^T$. The second replay has dot $-8$ and scalar $p=1$, so it sends the second coordinate to $-4-(-8)=4$. Every replay product, accumulated sum and subtraction is zero or one of the stipulated normal values, hence is represented exactly by [F5]. Thus the transformed RHS is $(1,4{)}^T$, and descending substitution with $R=-I_2$ has zero off-diagonal products and exact divisions by $-1$, giving $k_1=(-1,-4{)}^T$. The midpoint is $y+k_1\mathrm{/}2=(1\mathrm{/}2,0{)}^T$, so $F_1=J(y+k_1\mathrm{/}2)=(-1\mathrm{/}2,0{)}^T$; all these exact values are stored.

The second RHS is $F_1-k_1=(1\mathrm{/}2,4{)}^T$. The first stored reflector has dot $1$ and active $p=(1,0{)}^T$, so it sends the first coordinate to $1\mathrm{/}2-1=-1\mathrm{/}2$ and leaves the second at $4$. The second reflector has dot $8$ and scalar $p=1$, so it sends the second coordinate to $4-8=-4$. All replay products, sums and subtractions are again zero or stipulated normal represented values. Descending substitution with $R=-I_2$ is exact as in step 2.1 and gives $q_2=(1\mathrm{/}2,4{)}^T$. Hence $k_2=q_2+k_1=(-1\mathrm{/}2,0{)}^T$.

For either solve, $W=I_2$ makes $\widehat{v}=x$ and the computed residual zero, so the evaluator returns $\bar{r}=0$. The other norm upper bounds and ${\gamma }_2=2u\mathrm{/}(1-2u)<0.021$ show the first certificate expression is at most $0.021\cdot 2\cdot 5+0.01\cdot 10=0.31<1\mathrm{/}2$, so the evaluator can soundly return $\bar{W}=1\mathrm{/}2$ and then $\rho =1\mathrm{/}2$ for $\bar{r}+\bar{W}$. Thus both solves meet tolerance one. By [F1]–[F2], this successful call uses one factorization and two triangular solves.

Here $J=-100$ and $z=h(-100)=-100$. By [F1], $k=(-100)\mathrm{/}(1+100)=-100\mathrm{/}101$, hence $y_1=1-100\mathrm{/}101=1\mathrm{/}101$.

The backward-Euler equation [F4] is $z=1-100z$, so $101z=1$ and its unique scalar root is $z=1\mathrm{/}101$. This agrees with step 1.1; [F3] also identifies the Rosenbrock update as one Newton correction, which is exact here because the residual is affine.

Since $z=-100<0$, the factor in [F2] is $1\mathrm{/}101\in (0,1)$. For general decay parameter $\lambda >0$, the factor is $1\mathrm{/}(1+h\lambda )$, which tends to zero as $h\lambda \to \mathrm{\infty }$.

Here $W=1$ and $k_1=f(y)=-\lambda y$. Since $h{k}_1\mathrm{/}2=-2y$, the midpoint state is $-y$ and $F_1=f(-y)=\lambda y$.

The second right-hand side is $F_1-k_1=2\lambda y$, so $q_2=2\lambda y$ and $k_2=q_2+k_1=\lambda y$. Therefore $y_{+}=y+h\lambda y=5y$, giving amplification factor $5$.

For the exact-J branch, [F2] gives $R(-4)=[1-4(1-2\alpha )]\mathrm{/}(1+4\alpha {)}^2=(5-4\sqrt{2})\mathrm{/}(5-2\sqrt{2}{)}^2$, because $\alpha =1-1\mathrm{/}\sqrt{2}$. Since $25<32<36$, we have $5<4\sqrt{2}<6$, so $\mathrm{\mid }5-4\sqrt{2}\mathrm{\mid }<1$; also $2\sqrt{2}<4$, so $5-2\sqrt{2}>1$. Thus $\mathrm{\mid }R(-4)\mathrm{\mid }<1$. This differs from step 2.1 and proves that boundedness of $\tilde{J}$ alone does not preserve stiff damping.

The approximate-J calculation in step 2.1 and the exact-J factor comparison establish the stated damping failure. No order or convergence conclusion follows from this stiff-scaled calculation.