Exercises — Adaptive step-size — RK45, error control

(a) An embedded Runge–Kutta pair (RK45 / Fehlberg / Dormand–Prince). (b) $\epsilon =\Vert y^{(5)}-y^{(4)}\Vert$. For a system this is the tolerance-scaled RMS norm (accept when $\epsilon \le 1$); for a single equation it is just $|y^{(5)}-y^{(4)}|$. It is free because it reuses every expensive $f$-evaluation and adds only a cheap weighted difference $h{\sum }_i(b_i-b_i^{\ast })k_i$. (c) The controller uses $p=4$: the error estimate essentially measures the 4th-order result's local error. You still advance with $y^{(5)}$ (local extrapolation), but the exponent in the update is $1/(p+1)=1/5$ using $p=4$. See Order of a Numerical Method.

Rule: accept iff $\epsilon \le \text{tol}$.

• A: $2\times 10^{-7}\le 10^{-6}$ → accept.
• B: $10^{-6}\le 10^{-6}$ → accept (equality passes).
• C: $9\times 10^{-6}>10^{-6}$ → reject, redo the step from the same $t_n$.

(a) $\epsilon =5\times 10^{-8}<10^{-6}$ → accept. (b) Ratio $\frac{\text{tol}}{\epsilon }=\frac{10^{-6}}{5\times 10^{-8}}=20$. $h_{\text{new}}=0.9\times 0.20\times 20^{1/5}=0.18\times 1.8206=\mathrm{0.3277.}$ We had error to spare, so we grow — but only by $20^{1/5}\approx 1.82$, not by 20, because local error scales like $h^5$.

(a) $8\times 10^{-6}>10^{-6}$ → reject. (b) Ratio $=\frac{10^{-6}}{8\times 10^{-6}}=0.125$. $h_{\text{new}}=0.9\times 0.15\times 0.125^{1/5}=0.135\times 0.6598=\mathrm{0.08908.}$ Retry from the same $t_n$ with $h_{\text{new}}\approx 0.0891$.

(a) $\text{tol}=\text{atol}+\text{rtol}|y_n|=10^{-8}+10^{-6}\times 40=10^{-8}+4\times 10^{-5}\approx 4.001\times 10^{-5}$. (b) $\epsilon =3\times 10^{-5}<4.001\times 10^{-5}$ → accept (the relative part saved us). (c) Ratio $=4.001\times 10^{-5}/3\times 10^{-5}=1.3337$. $h_{\text{new}}=0.9\times 0.10\times 1.3337^{1/5}=0.09\times 1.0593=\mathrm{0.09534.}$ Note how using only atol ($10^{-8}$) would have rejected this perfectly good step — see Local vs Global Truncation Error for why absolute-only tolerances misbehave when $|y|$ is large.

Set $\epsilon =\text{tol}$. Then the raw update factor is $r=S{\left(\frac{\text{tol}}{\epsilon }\right)}^{1/5}=S\cdot 1=S=0.9,h_{\text{new}}=0.9h.$ Predicted error scales as $h^5$, so ${\epsilon }_{\text{next}}\approx C{h}_{\text{new}}^5=C(0.9h{)}^5=0.9^5(C{h}^5)=0.9^5\epsilon =0.9^5\text{tol}.$ $0.9^5=0.59049$, so ${\epsilon }_{\text{next}}\approx 0.590\text{tol}$ — comfortably under budget. That is exactly what the safety factor buys: even when we hit tol dead-on, the next step is predicted to land at ~59% of tol, so it is very likely accepted on the first try. This is the same trade-off drawn in the opening figure.

Shrink factor $=(\text{tol}/\epsilon {)}^{\text{exp}}=(1/4{)}^{\text{exp}}=0.25^{\text{exp}}$.

• Student X: $0.25^{1/5}=0.75786$ → $h$ drops to ~76% of old.
• Student Y: $0.25^{1/2}=0.5$ → $h$ halves. Because true local error $\propto h^5$, a factor-4 error overshoot only needs $h\to 4^{-1/5}h\approx 0.758h$ to fix. Student Y's $1/2$ exponent over-shrinks to $0.5h$, wasting steps: the retried step will come back at $\epsilon \approx C(0.5h{)}^5=4\cdot 0.5^5\text{tol}=0.125\text{tol}$ — 8× under budget, meaning the very next step immediately tries to grow again. Wrong exponent ⇒ oscillating, inefficient stepping. See Taylor Series Expansion for where the $h^5$ comes from.

(a) Ratio $=10^{-6}/2\times 10^{-13}=5\times 10^6$. $r_{\text{raw}}=0.9\times (5\times 10^6{)}^{1/5}=0.9\times 21.867=\mathrm{19.68.}$ (b) $19.68>f_{\max }=5$ → the clamp fires (recall $f_{\max }\approx 5$ from the opening definitions). (c) $h_{\text{new}}=h\cdot f_{\max }=0.04\times 5=0.20$. The $h^5$ error model is only trustworthy locally; a ~20× jump could leap into a region with completely different curvature, so we cap the growth at 5×.

Attempt 1: ${\epsilon }_1=6.25\times 10^{-6}>10^{-6}$ → reject. Stay at $t=0$. Ratio $=10^{-6}/6.25\times 10^{-6}=0.16$. Factor $=0.9\times 0.16^{1/5}=0.9\times 0.69314=0.62383$ (inside clamps). $h\to 0.10\times 0.62383=\mathrm{0.062383.}$ Attempt 2 (from $t=0$ with $h=0.062383$): ${\epsilon }_2=4\times 10^{-7}<10^{-6}$ → accept. Advance: $t=0+0.062383=0.062383$. Now pick the step for the next interval: ratio $=10^{-6}/4\times 10^{-7}=2.5$, factor $=0.9\times 2.5^{1/5}=0.9\times 1.2011=1.0810$. $h_{\text{next}}=0.062383\times 1.0810=\mathrm{0.067437.}$ Summary: rejected once (t stays 0), accepted second try (t advances to ≈0.06238), then grows slightly to ≈0.06744 for the next interval.

(a) Attempts $=80+12=92$. (b) $f$ calls $=92\times 6=552$. (c) Wasted $=12\times 6=72$ calls; fraction $=72/552=0.13043\dots \approx 13.0\%$. A rejection rate around 10–15% is healthy — it means the controller is pushing $h$ to the edge of the tolerance rather than being timidly small.

Attempt 1 ($h=0.20$): $\epsilon =3.2\times 10^{-5}>10^{-6}$ → reject. Ratio $=10^{-6}/3.2\times 10^{-5}=0.03125$. Factor $=0.9\times 0.03125^{1/5}=0.9\times 0.5=0.45$ (inside clamps, since $0.45>f_{\min }=0.2$). $h\to 0.20\times 0.45=\mathrm{0.09.}$ Attempt 2 ($h=0.09$): $\epsilon =3.5\times 10^{-6}>10^{-6}$ → reject. Ratio $=10^{-6}/3.5\times 10^{-6}=0.285714$. Factor $=0.9\times 0.285714^{1/5}=0.9\times 0.778087=0.700278$. $h\to 0.09\times 0.700278=\mathrm{0.063025.}$ Attempt 3 ($h=0.063025$): $\epsilon =5\times 10^{-7}<10^{-6}$ → accept. Advance: $t=0.063025$. Final: $t\approx 0.06303$ after two rejections and one acceptance. Notice the clamp never fired here; each shrink stayed well above $f_{\min }=0.2$ of the previous $h$.

(a) Ratio $=10^{-6}/0.3=3.3\overline{3}\times 10^{-6}$. Raw factor $=0.9\times (3.333\times 10^{-6}{)}^{1/5}$. $(3.333\times 10^{-6}{)}^{1/5}=0.080274$, so raw $=0.9\times 0.080274=0.072247$. (b) $0.072247<f_{\min }=0.2$ → the lower clamp fires. (c) $h_{\text{new}}=h\times f_{\min }=0.30\times 0.2=0.06$. (d) A raw factor far below $f_{\min }$ repeatedly is the fingerprint of a stiff region — an explicit adaptive method may thrash here; consider an implicit solver (see Stiff ODEs and Stability).

Ratio $=\text{tol}/\epsilon =32$.

• Wrong ($p=5$, exponent $1/6$): factor $=32^{1/6}=1.7818$, $h_{\text{new}}=0.17818$.
• Correct ($p=4$, exponent $1/5$): factor $=32^{1/5}=2.0$, $h_{\text{new}}=0.20$. The correct controller grows more (factor 2 vs 1.78). Reason: the error estimate $\epsilon$ tracks the 4th-order solution's error, which scales as $h^5$; a 32× error surplus permits a $32^{1/5}=2$× step increase. Using $p=5$ underestimates how fast you may grow — safe but wasteful. The mnemonic stays: "Tol Over Error, to the Fifth-root."

(a) Maximum shrink: $h_{\text{raw}}=h\cdot f_{\min }=8\times 10^{-6}\times 0.2=1.6\times 10^{-6}$. (b) $1.6\times 10^{-6}<h_{\min }=2\times 10^{-6}$ → it violates the floor. (c) The solver must not silently take a sub-floor step. Standard behaviour (e.g. scipy's solve_ivp) is to either clamp $h$ up to $h_{\min }$ and accept a slightly-too-large error with a warning, or abort with a "step size too small" failure. Blindly continuing to shrink is dangerous because (i) $h$ can underflow toward machine $\epsilon \approx 2.2\times 10^{-16}$, where $t_n+h$ rounds back to $t_n$ and the integration stalls forever; and (ii) tens of thousands of micro-steps accumulate round-off that swamps the very accuracy you were chasing. Repeatedly hitting $h_{\min }$ is the classic signal to switch to an implicit method built for stiff problems.