Exercises — Stiff equations — implicit methods, backward Euler

WHAT: just substitute. $z=h\lambda =(0.05)(-20)=-1$. WHY these numbers: $z$ bundles step size and decay into one dimensionless quantity — the only thing the stability test looks at. $R(z)=\frac{1}{1-z}=\frac{1}{1-(-1)}=\frac{1}{2}=\mathrm{0.5.}$ Since $|0.5|<1$, this step shrinks $y$: stable. ✅

$z=-1$, so $R(z)=1+(-1)=0$. $|0|<1$ → stable (in fact this hits the true value's fast decay hard). Recognising that forward Euler's factor is $1+z$ (not $1/(1-z)$) is the whole point of this rung.

Backward Euler. Its stable region (violet) is the entire outside of a small disk on the right, which contains the whole negative real axis — every negative $z$ is stable. Forward Euler's stable region (magenta) is a disk of radius $1$ centred at $-1$; it only covers $-2<z<0$. Beyond $z=-2$ forward Euler is unstable. This is the A-stability of backward Euler drawn as a picture.

WHAT: backward Euler on the linear test problem gives $y_{n+1}=\frac{1}{1-h\lambda }{y}_n$. WHY rearrange: $y_{n+1}$ appears inside $f$, so we move it to the left — the defining implicit move. $z=h\lambda =(0.5)(-4)=-2,R=\frac{1}{1-(-2)}=\frac{1}{3}.$ $y_1=\frac{1}{3}\cdot 1=0.3333,y_2=\frac{1}{3}\cdot y_1=\mathrm{0.1111.}$ Compare true $y(1)=e^{-4}=0.0183$ — decaying too, direction correct (first order, so not exact).

$z=h\lambda =(0.6)(-4)=-2.4$. Forward: $R=1+z=1-2.4=-1.4$. $y_1=-1.4$, $y_2=1.96$. $|-1.4|>1$ → grows and flips sign: unstable. Backward: $R=1/(1-(-2.4))=1/3.4=0.2941$. $y_1=0.2941$, $y_2=0.0865$ → decays smoothly: stable. The difference is entirely in the factor: forward left the disk of radius $1$; backward never can.

WHAT: $y_1=y_0+hf(t_1,y_1)=0+0.1(-30y_1+3)$. WHY: $f$ is evaluated at the new $y_1$, so $y_1$ sits on both sides — collect it. $y_1=-3y_1+0.3\Rightarrow y_1+3y_1=0.3\Rightarrow 4y_1=0.3\Rightarrow y_1=\mathrm{0.075.}$ (Sanity: the steady state is $y^{\prime }=0\Rightarrow y=3/30=0.1$; we moved from $0$ toward $0.1$. ✅)

Stability: $|1+h\lambda |<1\Rightarrow h<\frac{2}{|\lambda |}=\frac{2}{250}=0.008$. To reach $t=5$: $\text{steps}=\frac{5}{0.008}=625$ steps, just to stay stable, not for accuracy. Why this bound: the stable interval on the real axis is $-2<z<0$, i.e. $-2<h\lambda <0$, giving $h<2/|\lambda |$.

Let $z=h\lambda$. Since $h>0$ and $\lambda <0$, we have $z<0$, so $1-z=1-(\text{negative})=1+|z|>1$. Therefore $|R(z)|=\left|\frac{1}{1-z}\right|=\frac{1}{1-z}<1\text{because the denominator exceeds }1.$ No upper bound on $h$ appears anywhere — that is exactly what unconditional stability means. The step is chosen for accuracy alone.

Forward Euler must satisfy stability for the fastest mode: $h<2/|{\lambda }_2|=2/1000=0.002$. Steps $=10/0.002=5000$. Backward Euler has no stability cap; using $h=0.1$ gives steps $=10/0.1=100$. Ratio $=5000/100=50\times$ fewer steps. Why the fast mode dominates: stability of an explicit method is set by the smallest time scale $1/|{\lambda }_2|$, even long after that component is numerically zero — the definition of stiffness.

WHAT: the implicit equation is $g(y)=y-y_0-h(-y^3)=y-1+0.5y^3=0$. WHY Newton: simple iteration would need $|h\partial f/\partial y|<1$, the very stiff restriction we escape; Newton has no such cap. $g^{\prime }(y)=1+1.5y^2.$ At $y^{(0)}=1$: $g(1)=1-1+0.5=0.5$, $g^{\prime }(1)=1+1.5=2.5$. $y^{(1)}=y^{(0)}-\frac{g(y^{(0)})}{g^{\prime }(y^{(0)})}=1-\frac{0.5}{2.5}=1-0.2=\mathrm{0.8.}$ One more check of direction: $y$ should decay from $1$, and indeed $0.8<1$. ✅ (Further Newton steps refine toward the exact root of $g$.)

$z=h\lambda =(0.5)(-2)=-1$. Backward Euler: $y_1=\frac{1}{1-z}=\frac{1}{2}=0.5$. Error $=0.5-0.3679=0.1321$. Trapezoidal: factor $R=\frac{1+z/2}{1-z/2}=\frac{1-0.5}{1+0.5}=\frac{0.5}{1.5}=0.3333$. Error $=0.3333-0.3679=-0.0346$ — about $4\times$ smaller. Why: trapezoidal is second-order $O(h^2)$; backward Euler is first-order $O(h)$. Both are stable here, but trapezoidal is more accurate — a reminder that stability and accuracy are separate axes.

(a) Stiff? The stiffness ratio $=\frac{|{\lambda }_{\max }|}{|{\lambda }_{\min }|}=\frac{40000}{1}=40000$ — modes span 4 orders of magnitude of time scale. Yes, strongly stiff. (b) Forward Euler steps: governed by the fastest mode: $h<2/40000=5\times 10^{-5}$. Steps $=20/(5\times 10^{-5})=400000$. Prohibitive. (c) Choice: an L-stable implicit method such as backward Euler. Reason: the fastest mode has $z=h\lambda \to -\infty$ for any reasonable $h$; backward Euler's factor $1/(1-z)\to 0$ damps it in a single step, while an explicit method cannot. Choose $h$ purely for accuracy on the slow $\lambda =-1$ mode (e.g. $h\approx 0.05$–$0.1$), solving the implicit equations with Newton.

Define $E(h)=\frac{1/(1+h)-e^{-h}}{e^{-h}}=\frac{e^h}{1+h}-1$ (multiplying through by $e^h$). At $h=0.2$: $e^{0.2}=1.2214$, $1+h=1.2$, ratio $=1.01784$, so $E=0.01784\le 0.05$ ✅. At $h=0.3$: $e^{0.3}=1.34986$, $1+h=1.3$, ratio $=1.03835$, so $E=0.03835\le 0.05$ ✅. Both pass; the tolerance permits $h=0.3$ (and a little beyond). Why this matters: with stability free, the only thing sizing $h$ is this accuracy budget — exactly the situation the whole topic promised.

A-stability was shown in L3.2 (entire left half-plane stable). For the limit: ${\lim }_{|z|\to \infty }\left|\frac{1}{1-z}\right|={\lim }_{|z|\to \infty }\frac{1}{|1-z|}=0.$ The denominator grows without bound, so the factor vanishes. Both conditions hold → backward Euler is L-stable. Contrast with trapezoidal, whose factor tends to $|-1|=1$ (A-stable but not L-stable).

Everything reduces to one number ::: $z=h\lambda$: read $R(z)$, check $|R(z)|<1$, and remember backward Euler's $|R|<1$ for all $z$ in the left half-plane while forward Euler's holds only inside a disk of radius $1$ about $-1$.