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

Put $h=\Delta x$, let $C_h$ denote the sum of the four axial shifts, and let $K_h$ denote the sum of the four diagonal shifts. Reading both sides of the PDF stencil gives $M_h\dot{U}=L_{h}U$, where $M_h=2I/3+C_h/12$ and $L_h=(-10I/3+2C_h/3+K_h/6)/h^2$. The left-hand stencil acts on the time derivative; treating it as a scalar multiplier would lose the compact correction.

Write $D={\partial }_x^2+{\partial }_y^2$. Taylor expansion of the symmetric shifts givesandFor the exact diffusion solution $u_t=Du$, subtracting givesThe leading coefficient is not identically zero for general smooth data. Thus the normalized local residual is fourth order, and in the question's $h^{p+1}$ convention $\boxed{p=3}$. This is the compact semidiscrete nine-point diffusion scheme; its usual spatial numerical consistency order is four, rather than three.

Use the Fourier amplification symbol for the semidiscrete generator. At a mode $e^{i(k\xi +l\eta )}$, the two stencil symbols areThe mass symbol satisfies $1/3\le m\le 1$. The numerator of $\ell$ is increasing in each cosine, since each partial derivative is at least two, and its maximum is zero at $(1,1)$. Hence $\ell \le 0$ and the generator ${\lambda }_h=\ell /m$ is real and nonpositive.

Every exact-time mode is multiplied by $e^{t{\lambda }_h}$, of modulus at most one. Parseval identity therefore provesThe method is stable, uniformly in the spatial mesh. This is a semidiscrete conclusion: it imposes no time-step condition until a particular time integrator is applied. For the stated merely square-integrable initial data, cell averages give a well-defined stable initialization, with $\Vert U(0){\Vert }_{{\ell }_h^2}\le \Vert v{\Vert }_{L^2}$ by Cauchy–Schwarz. Arbitrary point sampling of an L2 equivalence class is not required.

For the linear multistep method, let $\rho (\zeta )={\zeta }^2-12\zeta /7+5/7$ and $\sigma (\zeta )=4{\zeta }^2/7-2/7$. The exponential-symbol order criterion for a multistep method compares $\rho (e^z)$ with $z\sigma (e^z)$. Direct expansion givesEquivalently, substitution of exact smooth values into the recurrence leaves $-2h^3{y}^{\prime \prime \prime }(t_n)/21+O(h^4)$, while the coefficients of $1,h,h^2$ cancel. Therefore $\boxed{\text{order }2}$; its unnormalized step defect is of order three. The nonzero cubic term excludes order three.

Factor the first characteristic polynomial:The unit root is simple and the other root is strictly inside the unit disk, so the method is zero-stable. Also $\rho (1)=0$ and ${\rho }^{\prime }(1)=\sigma (1)=2/7$, giving numerical consistency. The Dahlquist equivalence theorem therefore makes it convergent.

The mechanism can also be seen directly. The zero-stable homogeneous recurrence has bounded propagators because its two error modes are $1$ and $(5/7{)}^n$. Variation of constants for the perturbed recurrence bounds the error by a fixed multiple of its two starting errors and the sum of the step defects, plus $O(h)$ times a sum of errors from the Lipschitz vector field. Absorb the current implicit error term for sufficiently small $h$, then apply the discrete Gronwall inequality. On a fixed time interval the sum of $O(h^3)$ exact-start defects is $O(h^2)$, giving second-order numerical convergence with second-order starting values. The implicit step is locally well defined, for example by contraction when $4hL/7<1$ for a Lipschitz constant $L$.

For the Dahlquist test equation, write $z=h\lambda$. The amplification roots solveA root on $|\zeta |=1$ would require the boundary-locus valueIts denominator cannot vanish on that circle. Setting $\zeta =e^{i\theta }$ and simplifying the real part givesThus no amplification root can cross the unit circle while $\mathrm{Re}z<0$. The leading coefficient is also nonzero throughout that half-plane. At $z=-1$, the roots are $(6\pm \sqrt{3})/11$, both strictly inside the unit disk. Connect any other left-half-plane point to $-1$: continuity of the roots, without a boundary crossing or a root escaping through infinity, keeps both roots inside.

On the imaginary axis a boundary root can occur only at $\theta =0$, which gives $z=0$; there the roots $1$ and $5/7$ satisfy the simple-unit-root condition. Hence $\boxed{\text{the method is A-stable}}$, including the usual boundary root convention at zero.

Use the discrete inner product $⟨U,V{⟩}_h=h{\sum }_{m=1}^{J-1}{\overline{U}}_m{V}_m$, with endpoints $U_0=U_J=0$. The centered first-difference matrix $D_x$ is skew-Hermitian, while summation by parts givesSince $\alpha$ is real, the centered advection contribution has zero real part. ThereforeThis proves the energy contraction for centered drift-diffusion, and hence stability for every fixed real $\alpha$ and every mesh width, in the discrete L2 norm. The extra inequality $|\alpha |h\le 2$ would make the off-diagonal coefficients nonnegative and support a discrete maximum principle, but is not needed for this energy proof. No time-stepping method is involved here.

For a sufficiently smooth exact solution, the centered spatial residual isLet $e=U-R_{h}u$ be the nodal error. It satisfies $e^{\prime }=A_{h}e+{\tau }_h$, with $A_h=D_{xx}-\alpha D_x$ and $\Vert {\tau }_h{\Vert }_h\le C{h}^2$ on the fixed time interval. Part (a) gives $\Vert e^{t{A}_h}\Vert \le 1$, so Duhamel's formula yieldsThus the method is convergent, with second-order spatial error for smooth data and second-order compatible initialization. The homogeneous endpoint values introduce no boundary residual. For nonsmooth L2 data, stable approximating projections and a density argument extend numerical convergence, although the same order need not hold without the derivative bounds used above.

The nodes are $c=(0,1/2,1{)}^T$. Their Lagrange interpolation basis on the normalized step isIntegrating these polynomials from zero to each node gives $a_{ij}={\int }_0^{c_i}{\ell }_j(s)ds$. At $c_1=0$ the row vanishes; at $c_2=1/2$ the integrals are $(5/24,1/3,-1/24)$; at $c_3=1$ they are $(1/6,2/3,1/6)$. The endpoint integrals also give $b=(1/6,2/3,1/6{)}^T$. These are exactly the coefficients in the original PDF.

Indeed, the degree-three polynomial $q$ determined by $q(0)=y_n$ and $q^{\prime }(hs)={\sum }_j{\ell }_j(s)f(q(h{c}_j))$ satisfies $q^{\prime }(h{c}_i)=f(q(h{c}_i))$. Integrating it to the nodes gives precisely the stage equations, and integrating to the endpoint gives the update. This proves it is a collocation Runge-Kutta method, specifically the three-stage Lobatto IIIA method.

For an explicit order proof, the collocation integrals imply $A1=c$, $Ac=c^2/2$ and $A{c}^2=c^3/3$, with powers componentwise. The endpoint weights are Simpson weights, so $b^T{c}^j=1/(j+1)$ for $j=0,1,2,3$. Therefore the eight Butcher order conditions through order four are satisfied:For example the last expression is $b^T{c}^3/6$, because $A^{2}c=A{c}^2/2=c^3/6$. These conditions are obtained by matching all elementary differential terms through degree four in the exact and numerical Taylor expansions, not just by testing a scalar linear equation. Nonautonomous equations are included by adjoining $t^{\prime }=1$.

To exclude higher order, the scalar stability function computed below satisfies $R(z)-e^z=-z^5/720+O(z^6)$. The local defect is therefore nonzero at degree five even on the scalar test equation. Hence $\boxed{\text{the classical order is exactly }4}$.

Applying the stage equations to $y^{\prime }=\lambda y$ gives the stability functionThe poles are $3\pm i\sqrt{3}$, both in the right half-plane. To determine the whole domain rather than only its boundary, multiply numerator and denominator by 12 and put $z=x+iy$. Direct expansion givesThus $|R(z)|\le 1$ exactly when $x\le 0$, with no pole in that set. ConsequentlyOn the imaginary axis $|R|=1$, and in the open left half-plane it is less than one. Its limit along the negative real axis is one, so A-stability here does not imply stiff decay or L-stability.

Use $J$ equal intervals, $h=1/J$, and the interior piecewise-linear hat functions ${\varphi }_j$, which equal one at node $j$, vanish at all other nodes and vanish outside the two adjacent intervals. Write $u_h(x,t)={\sum }_{j=1}^{J-1}{U}_j(t){\varphi }_j(x)$. Testing with ${\varphi }_i$ and integrating the second derivative by parts givesThe integrals are local: $\int {\varphi }_i^2=2h/3$, $\int {\varphi }_i{\varphi }_{i+1}=h/6$, $\int ({\varphi }_i^{\prime }{)}^2=2/h$ and $\int {\varphi }_i^{\prime }{\varphi }_{i+1}^{\prime }=-1/h$. Nonadjacent functions do not overlap. Thus the consistent mass matrix and stiffness matrix areEquivalently the explicit nodal equations areFor merely L2 initial data a natural choice is the Galerkin projection $MU(0)=((v,{\varphi }_i){)}_i$; nodal interpolation is also possible for suitably smooth data. The relation $U^{T}MU=\Vert u_h{\Vert }_{L^2}^2$ identifies the natural finite-element energy. Differentiating gives $d(U^{T}MU)/dt=-2U^{T}KU\le 0$.

The sine vectors $v_j^{(m)}=\sin (j{\theta }_m)$, ${\theta }_m=m\pi /J$, simultaneously diagonalize these two symmetric tridiagonal matrices. Their eigenvalues areThus the semidiscrete decay rates are ${\lambda }_m=6(1-\cos {\theta }_m)/[h^2(2+\cos {\theta }_m)]$. The forward Euler method multiplies each mode by $1-k{\lambda }_m$. Stability requires $0\le k{\lambda }_m\le 2$ for every mode.

For a fixed finite Dirichlet grid the largest mode is ${\theta }_{J-1}=\pi -\pi /J$. Consequently the exact finite-grid bound isAs the mesh is refined, the upper bound decreases to $1/6$. Hence the mesh-independent answer is $\boxed{0\le \mu \le 1/6}$. This is the forward Euler limit for a consistent hat mass matrix. The frequently seen $1/2$ diffusion bound belongs to the standard second-difference scheme or a lumped mass matrix, and is not the bound for the consistent Galerkin equations derived here.

For constant-coefficient evolution on an infinite or periodic grid, spatial Fourier stability analysis reduces a multidimensional recurrence to independent modal problems. This gives exact L2 stability results when the coefficients and geometry admit that reduction; it also exposes high-frequency restrictions which a truncation-error calculation alone misses.

Consider a scalar one-step update $(S_{h,k}U{)}_j={\sum }_{\ell }{a}_{\ell }{U}_{j+\ell }$ on ${\mathbb{Z}}^d$. With $\hat{U}(\theta )={\sum }_j{U}_j{e}^{-ij\cdot \theta }$, its symbol is $G(\theta )={\sum }_{\ell }{a}_{\ell }{e}^{i\ell \cdot \theta }$ and ${\hat{U}}^{,n}=G^n{\hat{U}}^0$. Parseval identity givesThe upper bound is immediate by integration. For the reverse inequality choose Fourier data supported in a positive-measure neighborhood where $|G|$ approaches its essential supremum, then normalize them. This uses actual square-summable data, unlike an isolated infinite plane wave. On a periodic finite grid the discrete Fourier transform gives the same equality with a maximum over its resolved frequencies.

This proves the scalar Fourier power criterion. Contractivity is equivalent to $|G|\le 1$ at every frequency. For mesh-uniform stability on $nk\le T$, a sufficient condition is instead $\sup |G|\le 1+Ck$, because $(1+Ck{)}^n\le e^{CT}$. Conversely, if the powers at $n=\lfloor T/k\rfloor$ are bounded by a common $M_T\ge 1$, then $\sup |G|\le M_T^{1/n}=1+O(k)$ as $k\to 0$. Thus finite-time stability may allow mild modal growth; demanding strict contraction is stronger than necessary.

Stability connects numerical consistency to numerical convergence. Let $R_{h}u(t_n)$ be the exact solution restricted by a suitable projection and let its one-step defect be ${\tau }_n=R_{h}u(t_{n+1})-S_{h,k}{R}_{h}u(t_n)$. If $\Vert {\tau }_n\Vert \le Ck(h^p+k^q)$ and $\Vert S^n\Vert \le M_T$, the error recurrence givesThis proves numerical convergence for smooth solutions and convergent initial approximations. Density plus stability extends numerical convergence to the natural data space when stable projections and well-posed continuous evolution are available.

For completeness, the converse in the Lax equivalence theorem can be understood with faithful grid transfers. Let $J_h$ be a bounded reconstruction into the continuous data space and $R_h$ a bounded projection with $R_h{J}_h=I$, and suppose the grid and reconstructed norms are uniformly equivalent. numerical convergence for every datum, uniformly on the fixed time interval, makes the operators $J_h{S}_{h,k}^n{R}_h$ pointwise bounded. The Uniform boundedness principle makes their operator norms uniformly bounded. Apply this bound to $J_{h}V$, use $R_h{J}_{h}V=V$ and the norm equivalence, and obtain the uniform grid-power bound. With the numerical consistency estimate above, this proves the stability/numerical convergence equivalence under these standard linear assumptions.

Several examples illustrate the test. For forward time and centered space on the one-dimensional heat equation,Its minimum is $1-4\mu$, so $0\le \mu \le 1/2$ is the exact contraction range. In $d$ dimensions with mesh widths $h_a$, the symbol is $1-4{\sum }_a{\mu }_a{\sin }^2({\theta }_a/2)$, where ${\mu }_a=k/h_a^2$. Therefore ${\sum }_a{\mu }_a\le 1/2$, or $\mu \le 1/(2d)$ on an equal mesh. This follows by independently maximizing every sine factor, not by applying a one-dimensional bound separately to each direction.

For advection $u_t+a{u}_x=0$, forward time and centered space give $G=1-i\nu \sin \theta$, $\nu =ak/h$. At fixed nonzero advective Courant number its modulus exceeds one, so its powers grow without a mesh-uniform bound. Under the much smaller choice $k=O(h^2)$, however, $\log |G|\le a^2{k}^2/(2h^2)=O(k)$, allowing finite-time stability. This illustrates why the specified refinement path matters. For $a>0$, the upwind update instead hasand is contractive exactly for $0\le \nu \le 1$. Its numerical diffusion explains the stable damping as well as its first-order spatial accuracy.

For the heat equation with the Crank-Nicolson method, the symbol is $(1-2\mu {\sin }^2(\theta /2))/(1+2\mu {\sin }^2(\theta /2))$. It has modulus at most one for every nonnegative $\mu$, so there is no diffusion time-step restriction. Very stiff modes nevertheless have amplification near minus one, causing weakly damped oscillation; stability alone does not guarantee good stiff smoothing or small error.

For systems the symbol is a matrix, and the exact criterion is the uniform power bound for matrix Fourier symbols, not just a bound on eigenvalue moduli. For example $G=\left(\begin{array}{cc}1& 1\\ 0& 1\end{array}\right)$ has both eigenvalues one but $G^n=\left(\begin{array}{cc}1& n\\ 0& 1\end{array}\right)$ is unbounded as the number of steps increases. If $G=X\Lambda X^{-1}$ with uniformly bounded condition numbers, then $\Vert G^n\Vert \le \Vert X\Vert \Vert X^{-1}\Vert \max |\lambda {|}^n$, so a scalar-type condition does suffice. Multilevel schemes are treated by stacking their time levels into a companion matrix and checking its uniform powers. Simple unit-modulus roots and interior roots are the fixed-mode root condition, but uniformity near colliding roots still needs attention.

The technique is inexpensive, identifies the dangerous wavelengths, gives sharp constant-coefficient restrictions and generalizes directly to several space variables and Fourier-diagonalizable vector problems. Its limitations are equally concrete: boundaries can add unstable modes invisible to an infinite-grid calculation; variable coefficients couple frequencies; nonlinearities invalidate the fixed multiplier; nonuniform grids need another transform or energy estimate; and L2 stability does not by itself imply a maximum principle or positivity. Frozen-coefficient analysis is often useful evidence, while an energy or boundary estimate is needed to turn it into a global conclusion in those settings.

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.