Two-stage Radau IIA has order three and is scalar L-stable

The row sums of $A$ equal $c$, $b^T$ is the final row of $A$, and direct arithmetic gives ${\sum }_i{b}_i=1,b^{T}c=\frac{3}{4}\frac{1}{3}+\frac{1}{4}=\frac{1}{2},b^T{c}^2=\frac{3}{4}\frac{1}{9}+\frac{1}{4}=\frac{1}{3}\mathrm{.}$ Also $Ac=(1\mathrm{/}18,1\mathrm{/}2{)}^T$, so $b^{T}Ac=(3\mathrm{/}4)(1\mathrm{/}18)+(1\mathrm{/}4)(1\mathrm{/}2)=1\mathrm{/}6$.

Augment the equation by $Q(t,y)=(1,f(t,y))$. Because the rows of $A$ sum to $c$, the augmented stage times are exactly $t+c_{i}h$. Write $K_i=Q(X_i)$ for augmented stage derivatives. The invariant-ball bound [F3] gives $X_i-X=O(h)$, and the componentwise norm estimate [F4] yields, uniformly on the local tube, $K_i=Q+h{c}_i{Q}^{\mathrm{\prime }}Q+h^2({\sum }_j{a}_{ij}{c}_j{Q}^{\mathrm{\prime }}{Q}^{\mathrm{\prime }}Q+\frac{c_i^2}{2}{Q}^{\mathrm{\prime }\mathrm{\prime }}(Q,Q))+O(h^3)\mathrm{.}$ Indeed first $K_j=Q+O(h)$ gives $X_i-X=h{c}_{i}Q+O(h^2)$; substituting this into the first-order Taylor expansion gives the $h{c}_i{Q}^{\mathrm{\prime }}Q$ term, and substituting once more into the second-order Taylor expansion gives the displayed $h^2$ terms, with $O(h^3)$ remainder because $Q\in C^3$ and its third derivative is bounded.

On $y^{\mathrm{\prime }}=\lambda y$, direct solution of the two-stage linear system gives $\det (I-zA)=1-\frac{2z}{3}+\frac{z^2}{6},1+z{b}^T(I-zA{)}^{-1}\mathbf{1}=\frac{1+z\mathrm{/}3}{1-2z\mathrm{/}3+z^2\mathrm{/}6}\mathrm{.}$ Stiff accuracy makes the update the second stage, so this is the stated scalar stability function.

Multiplying the expansion of step 1.2 by $h{b}_i$ and summing gives $X+hQ+\frac{h^2}{2}{Q}^{\mathrm{\prime }}Q+\frac{h^3}{6}(Q^{\mathrm{\prime }\mathrm{\prime }}(Q,Q)+Q^{\mathrm{\prime }}{Q}^{\mathrm{\prime }}Q)+O(h^4),$ using all four identities in step 1.1. The exact solution has the same expansion: $X^{\mathrm{\prime }}=Q$, $X^{\mathrm{\prime }\mathrm{\prime }}=Q^{\mathrm{\prime }}Q$, and $X^{\mathrm{\prime }\mathrm{\prime }\mathrm{\prime }}=Q^{\mathrm{\prime }\mathrm{\prime }}(Q,Q)+Q^{\mathrm{\prime }}{Q}^{\mathrm{\prime }}Q$; applying scalar Taylor separately to each time-coordinate function, with the componentwise norm estimate in [F4] and the bounded fourth derivative obtained from $X^{\mathrm{\prime }}=Q(X)$ and $Q\in C^3$, gives remainder $O(h^4)$. Thus the local defect is $O(h^4)$.

Put $w=-z=s+it$, $s\ge 0$. Expanding real and imaginary parts gives ${\mid 1+\frac{2w}{3}+\frac{w^2}{6}\mid }^2-{\mid 1-\frac{w}{3}\mid }^2=2s+\frac{2}{3}{s}^2+\frac{2}{9}{s}^3+\frac{1}{36}{s}^4+(\frac{2s}{9}+\frac{s^2}{18})t^2+\frac{t^4}{36}\ge 0.$ If the denominator vanished, this inequality would force the numerator to vanish; the numerator $1-w\mathrm{/}3$ vanishes only at $w=3$, where the denominator is $1+2+9\mathrm{/}6\mathrm{\ne }0$. Hence there are no poles for $\mathrm{Re}z\le 0$, and the inequality gives $\mathrm{\mid }R(z)\mathrm{\mid }\le 1$ there. Since the numerator degree is one and denominator degree is two, $R(z)\to 0$ as $\mathrm{\mid }z\mathrm{\mid }\to \mathrm{\infty }$ in that half-plane. This is scalar L-stability.