Past exam of the mathematics course of the University of Cambridge / 2012 / iii / Paper 69 / Section II

For a discretized ordinary differential equation, error control must distinguish the local truncation error of one exact-start step from accumulated global error. For a smooth order-$p$ one-step method, these are usually $O(h^{p+1})$ and $O(h^p)$ respectively, provided perturbations propagate stably. A small local defect is useful only when the underlying problem and the numerical update do not amplify it excessively.

There are several practical local error estimators. In step-doubling error estimation, compute one step $Y_h$ and two half steps $Y_{h/2,h/2}$ from the same starting value. If the leading method error is $C{h}^{p+1}$, the two-half-step error is $C{h}^{p+1}/2^p$, up to higher terms. Their difference $D=Y_{h/2,h/2}-Y_h$ therefore givesfor exact minus fine-step error. Adding that estimate is Richardson extrapolation. This costs extra evaluations but gives an independent check on a single-method implementation; nonsmooth events can invalidate the assumed expansion.

An embedded Runge-Kutta pair shares stages and uses two output weight vectors to obtain orders $p$ and $q<p$. Their difference generally estimates the lower-order local error $O(h^{q+1})$, not the higher-order error directly. Stage reuse makes this much cheaper than independent integrations. Predictor-corrector error estimation uses the appropriately calibrated difference of two formulas; raw predictor/corrector differences require their known error constants. For multistep schemes, changing the step size also changes interpolation coefficients and requires controlled step ratios and proper history updates, rather than reusing constant-step coefficients.

Component scales should reflect the requested physical accuracy. An example of absolute and relative error tolerances isPositive absolute tolerances avoid singular scaling when a component vanishes, while relative tolerances track widely differing magnitudes. Accept the step if $E\le 1$; otherwise reject it without committing the new state or multistep history. If the estimated defect scales as $h^r$, a basic adaptive step-size controller proposesGrowth limits, a safety factor, and proportional-integral smoothing reduce erratic changes. The exponent must match the actual estimator: $r=p+1$ for a genuine order-$p$ local defect, but $r=p$ for the usual $p/(p-1)$ embedded difference. It is not always $1/p$ or always $1/(p+1)$. A zero estimate uses a finite maximum growth factor, not division by zero.

Local acceptance does not itself certify a global error tolerance. If a dense approximation $\eta (t)$ has residual $d(t)={\eta }^{\prime }(t)-f(t,\eta (t))$ and $f$ is Lipschitz continuous with constant $L$ on the relevant trajectory region, Gronwall's inequality gives the residual bound for global ODE errorThe discrete analogue propagates and sums local defects. A posteriori defect correction, repeat integrations with tighter tolerances, and adjoint-weighted estimates for a chosen final functional provide stronger global checks. A smaller local target may be needed in an unstable flow; trajectory sensitivity should not be mislabeled as a code defect.

For a stiff differential equation, absolute stability can constrain explicit steps far more tightly than accuracy alone. Use an appropriate implicit method and keep nonlinear stage/linear-solve errors below the discretization-error budget. A-stability, L-stability and B-stability address different linear damping and nonlinear contractivity properties. An error controller cannot repair a method outside its stability domain. Tolerance-driven adaptation also does not automatically preserve positivity, conservation laws or a symplectic structure; these may require a suitable method or an additional acceptance test.

Finally, event detection and dense output have their own errors. Locate an event with a consistent interpolant and account for its state sensitivity; a discontinuous right side may require restarting the integrator. Refining indefinitely eventually exposes roundoff error and conditioning limits. Reliable error control balances truncation, solution sensitivity, algebraic-solve error and floating-point error, while reporting failure to meet a target instead of endlessly shrinking a step.

Fourier stability analysis applies most directly to linear, constant-coefficient evolution schemes on uniform infinite or periodic grids. Translation invariance makes spatial shifts diagonal in discrete Fourier modes. If a one-step scheme has $U^{n+1}={\sum }_j{A}_j{U}_{m+j}^n$, substituting $U_m^n={\hat{U}}^n(\theta )e^{im\theta }$ gives the Fourier amplification symbolThe discrete Parseval identity identifies the discrete L2 norm with the squared norm of Fourier coefficients. Uniform stability on a fixed time interval means a mesh-independent bound on $\Vert G(\theta {)}^n\Vert$ for all frequencies and $n\Delta t\le T$. In the scalar case $|G|\le 1$ yields contraction. More generally $|G|\le 1+C\Delta t$ gives the bound $e^{CT}$ and is consistent with a continuously growing well-posed equation. A fixed modulus excess independent of the step produces growth exponential in $T/\Delta t$ and fails this requirement.

For systems, eigenvalue moduli alone are insufficient. The uniform power bound for matrix Fourier symbols also controls eigenvector conditioning and Jordan growth. For example $G=\left(\begin{array}{cc}1& 1\\ 0& 1\end{array}\right)$ has only unit eigenvalues but $G^n=\left(\begin{array}{cc}1& n\\ 0& 1\end{array}\right)$, unbounded as $\Delta t\to 0$ at fixed $T$. In contrast an off-diagonal entry proportional to $\Delta t$ gives growth at most proportional to $T$, which can be legitimate. A uniformly conditioned diagonalization is a sufficient reduction to scalar modal bounds; nonnormal matrices require direct power or energy estimates. Multilevel schemes similarly require the root condition for their frequency-dependent amplification polynomial and a uniform bound on the associated companion evolution, including limiting repeated-root cases.

The simplest parabolic example is the Forward Euler method with centered spatial diffusion. For $\mu =\Delta t/h^2$,is the sharp mesh-uniform contraction range. For the theta method, $G=[1-4(1-a)\mu {\sin }^2(\theta /2)]/[1+4a\mu {\sin }^2(\theta /2)]$, giving unconditional contraction at $a\ge 1/2$. Implicit stability permits larger steps, but accuracy still constrains them, and the Crank-Nicolson method does not strongly damp high-frequency rough data.

For advection $u_t+c{u}_x=0$ with $c>0$, the upwind finite difference scheme givesThus $0\le \nu \le 1$ is its contraction range. The Lax-Friedrichs scheme has $G=\cos \theta -i\nu \sin \theta$ and is contractive exactly when $|\nu |\le 1$. These express a Courant–Friedrichs–Lewy condition: the numerical domain of dependence must encompass the physical one for convergence. That necessary condition is not sufficient by itself. Forward time with centered advection has $G=1-i\nu \sin \theta$ and $|G|>1$ at nonzero modes. At a fixed nonzero Courant number it is unstable even if $|\nu |\le 1$. Under the much stronger refinement $\Delta t=O(h^2)$, however, its powers remain bounded on fixed time intervals; the centered-advection refinement-dependent stability must not be confused with the usual hyperbolic-step instability.

Symbols also quantify dispersion and dissipation. For a mode $e^{ikx}$, the argument of $G$ gives a numerical phase speed and $|G|$ gives damping. Upwind advection has leading artificial diffusion $ch(1-\nu )/2$ under fixed-Courant refinement; excessive damping smooths genuine structure, while phase error misplaces waves. Stability controls boundedness, not either error by itself. Smooth low-frequency Taylor expansions establish consistency; high-frequency modal bounds establish stability, so analyzing only fixed physical wavenumbers can miss grid-scale instabilities.

For a well-posed linear initial-value problem, the Lax equivalence theorem makes consistency plus a uniform stability estimate equivalent to convergence in the chosen norm. Boundaries and variable coefficients need extra care: Fourier diagonalization no longer handles the full update, and an interior or frozen-coefficient symbol does not prove boundary or global stability. Boundary energy estimates, a suitable matrix analysis or a separate boundary-stability argument must supplement it. Nonlinear schemes also require controls beyond this linear symbol test. These limitations delimit what the method proves rather than reducing its usefulness: in its proper setting, it turns an arbitrarily large grid evolution into exact and interpretable modal algebra.