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

Differentiate the squared Euclidean norm and use the skew-symmetric matrix identity $A^T=-A$:Thus $\Vert \mathbf{y}(t){\Vert }_2^2$ is constant, and continuity of the nonnegative square root gives

Set $\overline{\mathbf{y}}=({\mathbf{y}}_n+{\mathbf{y}}_{n+1})/2$. The implicit midpoint rule readsTaking the inner product with ${\mathbf{y}}_{n+1}+{\mathbf{y}}_n=2\overline{\mathbf{y}}$ givesThe last equality again follows from skew-symmetry. Mathematical induction therefore yields

Write $A_n=A({\mathbf{y}}_n)$. The trapezoidal rule isTaking the inner product with ${\mathbf{y}}_{n+1}+{\mathbf{y}}_n$ givesThe two quadratic terms vanish because each $A_j$ is skew-symmetric. Moreover ${\mathbf{y}}_{n+1}^T{A}_n{\mathbf{y}}_n=-{\mathbf{y}}_n^T{A}_n{\mathbf{y}}_{n+1}$. HenceUnlike the implicit midpoint rule, the trapezoidal rule evaluates $A$ at two different states, so the mixed terms need not cancel.

The characteristic polynomials of a linear multistep method areandFor the order conditions for a linear multistep method, put $\alpha =(a,-1-a,1)$ and $\beta =(-(1+5a)/12,,8(1-a)/12,,(5+a)/12)$. The defectsvanish for $q=0,1,2,3$, whileThereforeIndeed, at $a=-1$ one has $C_4=0$ but $C_5=-4/3\ne 0$.

Every member of the family has order at least three, so it is consistent. By the Dahlquist equivalence theorem, convergence is therefore equivalent to zero-stability, which is characterized by the root condition for a multistep method.

The roots of $\rho$ are $1$ and $a$. They lie in the closed unit disk exactly when $|a|\le 1$. At $a=1$, however, the unit root $1$ is repeated, whereas at $a=-1$ the distinct unit roots $1$ and $-1$ are both simple. Consequently

The Second Dahlquist barrier states that an A-stable linear multistep method has order at most two. Part a shows that every method in this family has order at least three, so no member can be A-stable.

The exceptional reducible case does not evade the conclusion. At $a=1$,Thus $\zeta =1$ is an amplification root for every $z=h\lambda$ with $\mathrm{Re}z<0$, whereas the multistep A-stability criterion requires all such roots to lie strictly inside the unit disk. Hence

Let $h=\Delta x$. A one-dimensional Taylor expansion givesAdding the three coordinate directions shows that the exact solution has stencil defectThis is the local defect of the seven-point Dirichlet Laplacian. Its inverse has max-norm size $O(h^{-2})$: this follows from the discrete maximum principle, for example by comparison with the grid restriction of $x_1(1-x_1)+x_2(1-x_2)+x_3(1-x_3)$, whose stencil is $-6h^2$. Therefore the grid error is $O(h^2)$.

The question writes this error as $O(h^{p+1})$, soIn the more usual terminology, the finite difference method is second-order accurate.

Each diagonal entry of $A$ is $-6$. Two distinct grid unknowns have matrix entry $1$ precisely when their nodes are nearest neighbours in one coordinate direction, and otherwise have entry $0$. Nearest-neighbour adjacency is symmetric: node $i$ is adjacent to node $j$ exactly when node $j$ is adjacent to node $i$. Hence $A_{ij}=A_{ji}$ for every pair, independently of the chosen ordering, and

Extend a grid vector $v$ by zero to the Dirichlet boundary. Pairing contributions along undirected nearest-neighbour edges gives the discrete energy identitywhere $E$ includes edges from an interior node to a boundary node. This is the three-dimensional version of summation by parts for the seven-point Dirichlet Laplacian.

The right side is nonpositive. If it vanishes, every pair of neighbouring values agrees; connectivity of the grid and the zero boundary values then imply $v=0$. Thus $v^{T}Av<0$ for every nonzero $v$, so $A$ is negative definite. In particular, zero is not an eigenvalue, and therefore

Take the real $L^2$ inner product of the equation with $u$. The homogeneous Dirichlet boundary conditions and integration by parts giveThus $\Vert u(t){\Vert }_{L^2}\le \Vert u_0{\Vert }_{L^2}$, and applying the same estimate to the difference of two solutions gives continuous dependence on the initial data. The constant-advection term is skew-symmetric under these boundary conditions and contributes no energy. Hence

Write the semidiscrete system as ${\mathbf{u}}^{\prime }=L_h\mathbf{u}+\alpha D_h\mathbf{u}$. With zero boundary values, the centered second-difference matrix $L_h$ is symmetric negative definite and the centered first-difference matrix $D_h$ is skew-symmetric. ThereforeThis is a mesh-uniform stability estimate for the centered convection-diffusion semidiscretization, valid for every real $\alpha$.

Both centered differences have local spatial error $O(h^2)$ for a sufficiently smooth solution. Stability plus consistency gives convergence, equivalently by the semidiscrete form of the Lax equivalence theorem. Thus

Let ${\Phi }_h$ be the one-step map of a numerical method. It is a time-symmetric numerical method when reversing the step exactly reverses the update:

Let ${\phi }_h$ be the exact flow map and suppose the method has order $p$ with a nonzero leading local truncation error:Inverting this expansion changes the sign of its leading perturbation, sowhere transport by the exact flow only changes $d$ by $O(h)$ and hence does not affect the leading parity. On the other hand, replacing $h$ by $-h$ in the first expansion givesTime symmetry equates these expressions, so $(-1{)}^{p+1}=-1$. Hence $p+1$ is odd and

The method has step map ${\Phi }_h(y)=y+hf(y)+h^{2}g(y)/2$. It is not time symmetric. It suffices to test the scalar linear equation $y^{\prime }=\lambda y$, for which the stability function isTime symmetry would require $R(-z)R(z)=1$, butfor generic $z$. Thus $\boxed{\text{method 1 is not time symmetric}}$.

The trapezoidal rule isInterchanging $y_n$ and $y_{n+1}$ while replacing $h$ by $-h$ leaves the equation unchanged. Its backward step is therefore exactly its inverse, and $\boxed{\text{method 2 is time symmetric}}$.

The Butcher tableau has one stage with $a_{11}=1/2$ and $b_1=1$, so it is the implicit midpoint ruleAgain, interchanging the endpoints and changing $h$ to $-h$ leaves the equation invariant. Hence $\boxed{\text{method 3 is time symmetric}}$.

A linear multistep method for $y^{\prime }=f(t,y)$ has the formIts characteristic polynomials of a linear multistep method are $\rho (\zeta )={\sum }_j{\alpha }_j{\zeta }^j$ and $\sigma (\zeta )={\sum }_j{\beta }_j{\zeta }^j$. Substitution of the exact solution, or equivalently expansion of $\rho (e^z)-z\sigma (e^z)$ at $z=0$, gives the order conditions for a linear multistep method. The method has order $p$ whenwith a nonzero coefficient at the next power. Its one-step defect is then $O(h^{p+1})$, while its global error is $O(h^p)$ when stability prevents accumulation from being amplified.

Order alone does not ensure convergence because the recurrence may contain growing parasitic modes. Zero-stability is the root condition for a multistep method: every root of $\rho$ lies in $|\zeta |\le 1$, and every root on the unit circle is simple. Consistency is the first-order conditionThe Dahlquist equivalence theorem states that a consistent linear multistep method is convergent if and only if it is zero-stable. For example, the leapfrog method $y_{n+2}-y_n=2h{f}_{n+1}$ has $\rho (\zeta )={\zeta }^2-1$; its roots $\pm 1$ are simple, so this second-order method is convergent, although its parasitic $-1$ mode can oscillate.

To study stiff decay, apply the method to $y^{\prime }=\lambda y$ and put $z=h\lambda$. The amplification factors are the roots ofThe linear stability domain consists of the $z$ for which these roots satisfy the strict interior condition appropriate away from the boundary, with simple unit roots where boundary stability is admitted. The method is A-stable when this domain contains the left half-plane, so every mode with $\mathrm{Re}\lambda <0$ remains bounded for every positive step size. The Second Dahlquist barrier says that an A-stable multistep method has order at most two.

The main families illustrate the tradeoffs. The explicit Adams-Bashforth methods interpolate past derivative values and are inexpensive, but their stability domains are bounded. The implicit Adams-Moulton methods include the new derivative value; the one-step second-order member is the trapezoidal rule, which is A-stable. The backward differentiation formulas interpolate past solution values and differentiate the interpolant; the first two are A-stable, while higher orders retain useful sectors of the left half-plane. Implicit methods require a nonlinear solve at each step, commonly initialized by an explicit predictor to form a predictor-corrector pair. Thus order controls consistency error, the root condition controls convergence as $h\to 0$, and the amplification polynomial controls whether a convergent method remains useful on stiff equations.

Assume $p\in L^{\infty }(-1,1)$ with $p(x)\ge p_0>0$, $q\in L^{\infty }(-1,1)$ with $q\ge 0$, and $f\in L^2(-1,1)$. The natural energy space is the Sobolev space $V=H_0^1(-1,1)$. Multiplyingby $v\in V$ and using integration by parts gives the weak formulationwhereThe form is bounded, and the Poincare inequality together with $p\ge p_0$ makes it coercive:The functional $\ell$ is bounded by the Cauchy-Schwarz inequality. The Lax-Milgram theorem therefore gives a unique weak solution.

Define the energyIf $y$ is the weak solution, then for every $w\in V$,Thus $y$ is the unique minimizer of $J$. Conversely, differentiating $J(y+tw)$ at $t=0$ recovers $a(y,w)=\ell (w)$, so the minimization and weak problems are equivalent.

The Ritz method chooses a finite-dimensional conforming space $V_h\subset V$ and minimizes $J$ over $V_h$. Its minimizer $y_h$ satisfiesFor a basis $\{{\varphi }_j\}$ this becomes $Kc=F$, withCoercivity makes the stiffness matrix $K$ symmetric positive definite, so the discrete minimizer exists and is unique.

Subtracting the continuous and discrete equations gives Galerkin orthogonality,In the energy norm $\Vert v{\Vert }_a=\sqrt{a(v,v)}$, for any $v_h\in V_h$,After cancellation and minimization over $v_h$, this proves the sharp symmetric form of the Céa lemma:

For a mesh $-1=x_0<x_1<\cdots <x_N=1$, take $V_h$ to be the continuous functions that are affine on each interval and vanish at the endpoints. The interior piecewise-linear hat functions form a basis, and each overlaps only its two neighbours, so $K$ is symmetric tridiagonal. In the illustrative case $p=1$ and $q=0$ on a uniform mesh, its diagonal and neighbouring entries areIf $y\in H^2(-1,1)$, its nodal interpolant $I_{h}y\in V_h$ satisfies $\Vert y-I_{h}y{\Vert }_{H^1}=O(h)$. The best-approximation estimate therefore gives $\Vert y-y_h{\Vert }_{H^1}=O(h)$, and a standard duality argument improves the $L^2$ error to $O(h^2)$ under the corresponding elliptic regularity. This completes the route from the boundary-value problem through its variational principle to a sparse, stable finite-element approximation.