Past exam of the mathematics course of the University of Cambridge / 2012 / iii / Paper 69 (Numerical Solution of Differential Equations)

Let $I=(-1,1)$ and use the clamped second-order Sobolev spaceHere weak derivatives up to order two are square integrable, with Sobolev norm $\Vert v{\Vert }_{H^2}^2=\Vert v{\Vert }_2^2+\Vert v^{\prime }{\Vert }_2^2+\Vert v^{\prime \prime }{\Vert }_2^2$. The boundary values are the continuous traces given by the Sobolev trace theorem, not arbitrary pointwise representatives. This space is the $H^2$ closure of $C_c^{\infty }(I)$, so it is a closed Hilbert space. The derivative boundary conditions are essential, and $H_0^1(I)$ alone is insufficient.

Define the symmetric bilinear form and bounded linear functionalThen the variational problem isFor complex data replace $a$ by its sesquilinear version, take the real part of the forcing term and use absolute squares in $J$; the same argument applies.

Its first variation is $J^{\prime }(u)[v]=a(u,v)-\ell (v)$. Thus stationary points satisfy $a(u,v)=\ell (v)$ for every $v\in H$, the weak formulation of a clamped fourth-order equation. For smooth coefficients and functions, twice integrating the highest term and once the middle term gives the differential equation, with all boundary terms zero because $v=v^{\prime }=0$ at both ends. With only the stated continuous coefficients, this manipulation is instead understood distributionally: compactly supported test functions give $(p{u}^{\prime \prime }{)}^{\prime \prime }-(q{u}^{\prime }{)}^{\prime }+ru=f$ as a distribution. This avoids assuming nonexistent classical derivatives of $p$ or $q$.

By the extreme value theorem, $p_0={\min }_{[-1,1]}p$ is strictly positive, and $p,q,r$ are bounded. The Cauchy-Schwarz inequality therefore givesFor completeness, the clamped interval derivative norm bound follows from $v^{\prime }(-1)=v(-1)=0$ and the one-dimensional Sobolev representative: integrate $v^{\prime \prime }$ to recover $v^{\prime }$, and then $v^{\prime }$ to recover $v$. Applying Cauchy-Schwarz and integrating over an interval of length two givesConsequently $\Vert v{\Vert }_{H^2}^2\le 21\Vert v^{\prime \prime }{\Vert }_2^2$. Nonnegativity of $q,r$ yields the precise coercive bilinear form boundThe Lax-Milgram theorem states that a bounded bilinear form on a real Hilbert space satisfying $a(v,v)\ge c\Vert v{\Vert }^2$ with $c>0$ represents every bounded linear functional uniquely: there is exactly one $u$ with $a(u,v)=\ell (v)$ for all $v$. All of its hypotheses have now been verified on $H$.

This weak solution is the unique global minimum, rather than merely a stationary point. The Ritz-Galerkin equivalence for a symmetric coercive form here is the exact identityIt is strictly positive for every nonzero $w\in H$. Conversely, a minimum has zero first variation in every direction, so it satisfies the same weak equation. Density of compactly supported test functions and boundedness of $a,\ell$ extend the distributional equation back to all $H$ test vectors. There is exactly one minimum, and it is exactly the clamped weak solution. No additional coefficient smoothness or classical fourth derivative is needed.

The coefficients are those of the three-stage Lobatto IIIA method. Set $e=(1,1,1{)}^T$, $c=Ae=(0,1/2,1{)}^T$ and $C=\mathrm{diag}(c)$. Direct calculation gives all eight fourth-order conditions for a Runge-Kutta method:Powers of $c$ are componentwise. These Butcher order conditions prove order at least four for a smooth general nonlinear vector field; matching a scalar test equation alone would not suffice. A necessary fifth-order condition fails:Hence the method has order exactly four, with one-step local truncation error $O(h^5)$ and, under stability and smoothness hypotheses, global error $O(h^4)$.

For the Dahlquist test equation with $z=h\lambda$, the stage equations give the stability function $R(z)=1+z{b}^T(I-zA{)}^{-1}e$. Substitution simplifies it toThe denominator zeros are $3\pm i\sqrt{3}$, so no pole lies in the closed left half-plane. Moreover,Thus $|R(z)|\le 1$ throughout that half-plane: the method is A-stable. On the imaginary axis its amplification modulus is one. Since $R(z)\to 1$ as $z\to -\infty$, it is not L-stable; linear stability does not guarantee strong damping of very stiff modes.

Algebraic stability of a Runge-Kutta method requires $b_i\ge 0$ and a positive-semidefinite symmetric matrixAlthough all weights are positive, its first diagonal entry isThe quadratic form on the first coordinate vector is negative, so $M$ cannot be positive semidefinite. The method is not algebraically stable. This distinguishes nonlinear contractivity from the A-stability proved in part 2.

Let $h=\Delta x$ and sample a sufficiently smooth periodic exact solution. The central finite differences have expansionsandThe semidiscrete right side therefore equals $u_{xx}+\alpha u_x+h^2(u_{xxxx}/12+\alpha u_{xxx}/6)+O(h^4)$. The spatial consistency of a numerical method is second order, for each fixed $\alpha$. There is no time discretization in the stated ordinary differential system, so no separate temporal approximation order is being assigned. The Taylor statement concerns smooth solutions; L2-compatible initialization of grid data handles nonsmooth initial data without undefined point sampling.

Use the discrete norm $\Vert U{\Vert }_h^2=h{\sum }_{m=1}^N|U_m{|}^2$, with periodic indices. Two independent arguments prove uniform stability of a numerical method.

For Fourier stability analysis, a discrete Fourier mode $U_m(t)={\hat{U}}_k(t)e^{im{\theta }_k}$ with ${\theta }_k=2\pi k/N$ has semidiscrete eigenvalueIts real part is nonpositive, so $|{\hat{U}}_k(t)|\le |{\hat{U}}_k(0)|$. The unitary discrete Fourier transform and discrete Parseval identity then give $\Vert U(t){\Vert }_h\le \Vert U(0){\Vert }_h$ for every $N$, $t\ge 0$ and fixed real $\alpha$.

For the energy method, denote the centered first difference by $D_0$ and the forward difference by $D_{+}$. Periodic discrete summation by parts makes $D_0$ skew-adjoint and gives $⟨U,D_{xx}U{⟩}_h=-\Vert D_{+}U{\Vert }_h^2$. ThereforeThe constant mode is preserved rather than strictly damped. Both proofs apply equally to differences of two solutions. The semidiscretization is contractive in the discrete L2 norm. A later choice of a time integrator would require its own stability analysis; no time-step condition is needed for this continuous-time result.

The characteristic polynomials of a linear multistep method areExpanding the exponential-symbol order criterion for a multistep method givesThus the formal local consistency order is two for $\alpha \ne 2$, and three at $\alpha =2$, since the fourth-degree coefficient there is $-1/12$. This exceptional formal order is degenerate: ${\rho }^{\prime }(1)=\sigma (1)=0$, so it is not a convergent third-order time integrator.

For zero-stability, factor $\rho (\zeta )=(\zeta -1)(\zeta -(\alpha -1))$. The root condition for a multistep method requires every root to have modulus at most one and every unit-modulus root to be simple. For real $\alpha$ it holds exactly forAt zero the second root is $-1$, simple and admissible; at two the root one is double and inadmissible. Outside the closed interval the second root has modulus greater than one. The Dahlquist equivalence theorem says that a consistent linear multistep method is convergent exactly when it is zero-stable, for a locally Lipschitz initial-value problem with consistent starting values and solvable implicit steps. Hence the displayed interval is also exactly the convergence range, and every convergent member has order two.

The common-factor cancellation defect in a multistep recurrence explains why the $\alpha =2$ exception cannot be repaired just by canceling $\zeta -1$. For $y^{\prime }=0$, its recurrence allows linear growth from a vanishing but insufficiently small starting perturbation: $y_0=0$, $y_1=\sqrt{h}$ gives $y_n=n\sqrt{h}$, diverging at a fixed positive time. Canceling the factor silently imposes an additional starting relation and changes the method.

Every convergent member is also A-stable. For $0<\alpha <2$, put $\zeta =e^{i\theta }$ and calculate its boundary-locus test for multistep A-stability:The denominator never vanishes on the unit circle in this parameter interval. It follows thatTherefore no amplification root can cross the unit circle as $z=h\lambda$ moves in the open left half-plane. For small negative real $z$, the root near one moves to $1+z+O(z^2)$, inside the circle; the other begins at $\alpha -1$, also strictly inside. The leading coefficient $1-(1-\alpha /4)z$ never vanishes in the left half-plane. Continuity of polynomial roots and connectedness of that half-plane now keep both roots inside everywhere. On its imaginary boundary the same locus has zero real part only at $\zeta =1$, corresponding to $z=0$; that root is simple because $2-\alpha >0$.

At $\alpha =0$, the two interlaced subsequences follow the trapezoidal rule with step $2h$:The right side has modulus at most one for $\mathrm{Re}z\le 0$. Its unit-circle roots are simple; the double zero at $z=-1$ lies strictly inside and is harmless. Thus the endpoint is also A-stable. Combining this with part 1, the requested range is $0\le \alpha <2$, including the BDF2 method at $\alpha =4/3$.

Take the endpoint values $u_0=u_M=0$ and evolve the interior indices $1,\dots ,M-1$. This matches $\Delta x=1/M$ on $[0,1]$. If the printed upper index $m=M$ is retained, the odd ghost extension $u_{M+1}=-u_{M-1}$ makes its boundary equation identically compatible; it does not add an interior unknown.

Writing $h=\Delta x$ and $\tau =\Delta t$, this is the theta method with parameter $a$:For a smooth exact heat solution, its normalized residual at $t_n$ isThus the generic consistency orders are first in time and second in space for $a\ne 1/2$, and second in both time and space for $a=1/2$, the Crank-Nicolson method. These are independent-step orders. Along the special parabolic refinement $\tau =\mu h^2$, the leading combined residual cancels if $a=1/2-1/(12\mu )$, giving a higher combined smooth-data order; it does not change the generic temporal order at independently varied meshes.

The given L2 initial condition need not have pointwise values or four derivatives. Initialize by a bounded cell-average or projection operator, as in L2-compatible initialization of grid data. The Taylor estimate applies to smooth data or sufficiently positive times after heat smoothing, rather than an unqualified uniform high-order estimate at the initial instant for arbitrary rough data. Stability and convergence for rough data are distinct from this consistency order.

The Dirichlet difference Laplacian has an orthogonal sine basis with eigenvalues $-{\eta }_j/h^2$, where ${\eta }_j=4{\sin }^2(j\pi /(2M))$ for $j=1,\dots ,M-1$. On a mode the theta method amplification isIf $a\ge 1/2$, the denominator is positive andThus $|G(s)|\le 1$ for every mode and every $\mu >0$. Orthogonality gives a uniform discrete L2 norm contraction, not merely individual eigenvalue estimates.

If $a<0$, a positive $s=-1/a$ makes the implicit matrix singular. If $a=0$, $G=1-s$ becomes arbitrarily large. If $0<a<1/2$, $G(s)$ tends to $-(1-a)/a$, whose magnitude is greater than one. Hence unconditional stability holds exactly for $a\ge 1/2$.

For completeness, $a<1/2$ is contractive under $\mu {\eta }_{\max }\le 2/(1-2a)$; the mesh-independent sufficient restriction is $\mu \le 1/[2(1-2a)]$. This follows from the same squared-modulus difference and keeps the denominator positive even for negative $a$. At $a=1/2$, very high frequencies approach amplification $-1$, so A-stability does not imply damping of rough heat data. At $a=1$ the Backward Euler method damps them strongly; this is the rough-data damping distinction for the theta method.

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.