Runge–Kutta Methods and Order — Examples

On the first midpoint stage, $k_1=f(-h\mathrm{/}2,0)=0$, so the second-stage state remains $0$. Its time is $-h\mathrm{/}2+h\mathrm{/}2=0$, where the definition of $f$ gives $k_2=f(0,0)=0$. Thus the method output is $0$.

The displayed piecewise function satisfies $y(0)=0$ and $y^{\mathrm{\prime }}(t)=f(t,y(t))$ for $t<0$ and $t>0$. At $t=0$ its two one-sided derivatives are both $0=f(0,y(0))$, so it is an exact solution there as well. Its value at $h\mathrm{/}2$ is $h^2\mathrm{/}8$, hence the local defect is $0-h^2\mathrm{/}8=-h^2\mathrm{/}8$.

If one constant $C$ gave the uniform order-two estimate in [F3], these admitted exact-start steps would satisfy $h^2\mathrm{/}8\le C{h}^3$, hence $1\mathrm{/}(8h)\le C$ for every $h>0$. The left side is unbounded as $h\downarrow 0$, so no such uniform constant exists. The field is continuous but not differentiable in time at $0$, identifying the missing smoothness condition.

By [F4]–[F5], $y(t)=\exp (-t)$ solves $y^{\mathrm{\prime }}=-y$ with $y(0)=1$. For one step, the tableaus [F1] and factors [F2] give $R$; put $d_h:=R-e^{-h}$. Taylor's theorem [F3] applied to $e^{-h}$ gives, for $p=1,2,4$ in the respective methods, $d_h=c_0{h}^{p+1}+O(h^{p+2}),c_0=-\frac{1}{2},\frac{1}{6},\frac{1}{120}\mathrm{.}$ The first omitted exponential terms after the Euler, midpoint, and RK4 polynomials in [F2] are respectively $+h^2\mathrm{/}2$, $-h^3\mathrm{/}6$, and $-h^5\mathrm{/}120$; subtracting them from those polynomials gives the listed defect coefficients.

For each method and sufficiently small $h$, set $q_h:=R{e}^h=1+e^h{d}_h$, using $e^{-h}{e}^h=e^0=1$ from [F6]. Thus $\mathrm{\mid }{q}_h-1\mathrm{\mid }\le C{h}^{p+1}$ for a constant independent of $N$. For $m\le N=1\mathrm{/}h$, factoring $q_h^m-1$ gives $\mathrm{\mid }{q}_h^m-1\mathrm{\mid }\le m\mathrm{\mid }{q}_h-1\mathrm{\mid }\max (1,\mathrm{\mid }{q}_h\mathrm{\mid }{)}^{m-1}\le C{h}^p{e}^C\mathrm{.}$ The last bound follows from $\mathrm{\mid }{q}_h\mathrm{\mid }\le 1+C{h}^{p+1}\le e^{C{h}^{p+1}}$ by [F5], so the factor of powers is at most $e^{CN{h}^{p+1}}\le e^C$.

Repeated use of [F6] gives $(e^{-h}{)}^N=e^{-Nh}=e^{-1}$ and $e^{-h(N-1)}=e^{-1+h}$. Since $R=q_h{e}^{-h}$, the finite-factor identity yields $R^N-e^{-1}=(R-e^{-h}){\sum }_{k=0}^{N-1}{R}^{N-1-k}{e}^{-kh}=d_h{e}^{-h(N-1)}{\sum }_{m=0}^{N-1}{q}_h^m\mathrm{.}$ By step 2.1 each summand is $1+O(h^p)$ uniformly in $m$, and $e^{-1+h}=e^{-1}(1+O(h))$ by [F3]–[F4]. Thus the sum factor is $N{e}^{-1}(1+O(h))$.

Substituting the expansion from step 1.1 and $N=1\mathrm{/}h$ into step 3.1 gives $R^N-e^{-1}=c_0{e}^{-1}{h}^p(1+O(h))\mathrm{.}$ Because $c_0\mathrm{\ne }0$, taking absolute values gives the three positive constants and slopes in the Example.

The eigenvalues of $A$ are $-1$ and $-200$, so $\rho =200$ and $h\rho =100$. For $s=8$, the certified RKC1 bound is $2s^2=128$, hence the full step is stable on both eigenmodes and the matrix update is nonexpansive.

The corresponding RKC1 scalar arguments are $z=-1\mathrm{/}2$ and $z=-100$; the latter maps to $1+z\mathrm{/}{s}^2=-9\mathrm{/}16\in [-1,1]$. The first-order recurrence uses eight stage evaluations for this step.

For Euler on the $-200$ mode, one step of length $1\mathrm{/}2$ has $z=-100$ and factor $1+z=-99$, whose absolute value exceeds one. If Euler instead divides the same time interval into $m$ equal substeps, stability on that mode requires $(h\mathrm{/}m)200\le 2$, hence $m\ge 50$. For $m\ge 50$, the Euler argument on the $-1$ mode has magnitude $h\mathrm{/}m=1\mathrm{/}(2m)\le 1\mathrm{/}100\le 2$, so both modes are stable. Thus the eight-stage RKC1 step uses fewer right-hand-side evaluations than stable Euler substepping for this spectrum, without implying comparable accuracy.

The matrix is nonnormal: $A{A}^{\mathsf{T}}=\left(\begin{array}{cc}17& -4\\ -4& 1\end{array}\right)$ while $A^{\mathsf{T}}A=\left(\begin{array}{cc}1& -4\\ -4& 17\end{array}\right)$. Its two eigenvalues are $-1$, so those of $hA$ are both $-1\mathrm{/}2\in [-2,0]$ and the scalar stability factor is $1\mathrm{/}2$.

The one-stage update is $I+hA=\left(\begin{array}{cc}1\mathrm{/}2& 2\\ 0& 1\mathrm{/}2\end{array}\right)\mathrm{.}$ It sends $v=(0,1)$ to $(2,1\mathrm{/}2)$, whose Euclidean norm is $\sqrt{17}\mathrm{/}2>1=\mathrm{\parallel }v{\mathrm{\parallel }}_2$. Therefore eigenvalue-wise scalar stability does not imply Euclidean nonexpansion here. However, for $u=(1,1)$, $⟨Au,u⟩=2>0$, so this matrix also violates the nonpositivity hypothesis in [F2].

To test omission of self-adjointness while retaining nonpositivity, take $B=\left(\begin{array}{cc}-1& 2\\ 0& -1\end{array}\right)$ and $\tilde{h}=1$, still with $s=1$. For every $u=(x,y)$, $⟨Bu,u⟩=-x^2+2xy-y^2=-(x-y{)}^2\le 0$. Both eigenvalues of $B$ are $-1$, so $\rho =1$, $\tilde{h}\rho =1\le 2$, and $\tilde{h}{\lambda }_i=-1\in [-2,0]$. Yet $B^{\mathsf{T}}\mathrm{\ne }B$ and $I+\tilde{h}B=\left(\begin{array}{cc}0& 2\\ 0& 0\end{array}\right)$ sends $(0,1)$ to $(2,0)$, increasing its Euclidean norm from $1$ to $2$. Thus the self-adjoint hypothesis in [F2] cannot be omitted even when nonpositivity and the spectral step restriction are retained.

Since $A=0$, the unperturbed stages all equal $y$. If $e_j$ is the perturbed-minus-unperturbed stage difference, then $e_0=0$, $e_1=ϵv$, and [F1] gives $e_j=2e_{j-1}-e_{j-2}+ϵv$ for $2\le j\le s$.

Induction gives $e_j=ϵj(j+1)v\mathrm{/}2$: it holds for $j=0,1$, and substituting the formulas at $j-1,j-2$ into the recurrence yields $ϵ\{(j-1)j-(j-2)(j-1)\mathrm{/}2+1\}v=ϵj(j+1)v\mathrm{/}2$.

At the output $j=s$, $\mathrm{\parallel }{e}_s{\mathrm{\parallel }}_2=ϵs(s+1)\mathrm{/}2$ because $\mathrm{\parallel }v{\mathrm{\parallel }}_2=1$. The assumptions of [F2] hold for $A=0$, so the lemma supplies the matching upper bound; this perturbation pattern attains equality in its quadratic stage-count estimate.