Past exam of the mathematics course of the University of Cambridge / 2005 / iii / Paper 69 / 7

An ODE solver must distinguish the error made by one step started from exact data from the error accumulated along the computed trajectory. For an order-$p$ one-step map ${\Phi }_h$, define the exact-start defect $d_n=y(t_{n+1})-{\Phi }_{h_n}(t_n,y(t_n))=O(h_n^{p+1})$. If the step map is Lipschitz with factor $1+L{h}_n$, its global error obeysIteration, using $1+Lh\le e^{Lh}$, proves the variable-step global ODE error boundIt explains the difference between local defect order $p+1$ and global order $p$, and why local tolerance control alone does not certify a prescribed global error. Strong dynamical amplification can magnify many individually small defects.

A practical solver needs an inexpensive local error estimator. In step-doubling error estimation, compare one full step $Y_h$ with two half steps $Y_f$. Smoothness and propagation through the second half step give $Y_h=y_{\mathrm{exact}}+C{h}^{p+1}+O(h^{p+2})$ and $Y_f=y_{\mathrm{exact}}+C{h}^{p+1}/2^p+O(h^{p+2})$. ThereforeThe difference estimates the finer local correction with the correct denominator and sign; adding it cancels the leading error by Richardson extrapolation. The common starting value and the smooth error expansion are essential, and the extra evaluations make this estimator relatively expensive.

An embedded Runge-Kutta pair instead shares the stage computations and uses two output weight vectors. For orders $p$ and $p-1$, their difference is generally $O(h^p)$ and estimates the lower-order local defect, even when the higher-order output is accepted. Its controller exponent is therefore $1/p$, not $1/(p+1)$. For example the explicit Euler and Heun outputs share the first slope; their difference is $h(k_2-k_1)/2=O(h^2)$, while Heun's own defect is $O(h^3)$. Predictor-corrector differences can also estimate a defect when scaled by the appropriate leading error constants; an uncalibrated raw difference need not be the corrected method's error.

For a vector equation, choose positive component scalesThe absolute and relative error tolerances allow meaningful tests both near zero and for large components with different physical scales. The max-norm criterion $E\le 1$ controls every scaled estimated component. A root-mean-square criterion is another common choice but has a weaker individual-component implication.

If the estimate scales like $E\simeq C{h}^q$, an adaptive step-size controller useswith lower and upper bounds on the change factor. A rejected step leaves the old state and time unchanged and recomputes with a smaller step; an accepted step advances with the chosen output and adjusts the next step. A zero estimate receives capped finite growth rather than division by zero. The safety factor allows for variation in the error constant. Step doubling has $q=p+1$, whereas a $p/(p-1)$ embedded pair usually has $q=p$. If a solver instead controls defect per unit time, division by $h$ changes the scaling exponent and must be included in the controller design.

Using previous estimates can smooth the sequence of accepted steps. A proportional-integral controller can take $h_{n+1}/h_n=\eta E_n^{-a}{E}_{n-1}^b$ with $0<b<a$, appropriate scaling by the estimator order and bounded factors. It avoids overreacting to a single noisy estimate, but needs safeguards after rejected steps, changes of order and abrupt coefficient changes. Event locations and discontinuities require locating the event and restarting the error expansion; blindly applying a smooth-step estimator across a jump can be misleading.

For a multistep method, adaptive step sizes change the interpolation history and coefficients. Constant-step coefficients cannot simply be reused with arbitrary step ratios. Variable-step formulas, compatible starting/restarting procedures and bounded ratios preserve the needed stability; local error and sometimes order are selected together to balance evaluation cost against accuracy.

Stiffness imposes an additional requirement. On a rapidly decaying mode $y^{\prime }=-\kappa y$, explicit Euler requires $h\kappa \le 2$ for absolute stability. Once $y$ is small, a local absolute-tolerance test can propose a large step even though roundoff or a small perturbation would then be amplified. A stiff solver therefore also uses a suitable implicit A-stable or L-stable method and solves its algebraic equations accurately enough for the estimator to be meaningful. A-stability alone does not ensure stiff decay, as the Lobatto stability function in question 4 tends to one for large negative arguments.

Finally, over-tightening tolerance eventually encounters roundoff and algebraic-solve error rather than further truncation-error improvement. Inexact Newton or linear solves must have errors below the local budget; a tiny time step with an inaccurate stage solve does not deliver the formal order. Global error assessment can use controlled refinements, defect propagation or adjoint weighting for a chosen output. Sensible error and step-size control therefore combines asymptotic estimators, tolerance scaling, stability, event handling and implementation error, rather than treating the accepted local estimate as a complete guarantee about the trajectory.