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

Expand about $t=t_{n+2}$ and use $y^{\prime \prime }=f^{\prime }(y)f(y)$ from the chain rule. The residual of the exact solution isThe coefficients through $h^3$ vanish, and the displayed fourth-order coefficient does not. Thus the multiderivative multistep method has order three.

At $h=0$, its first characteristic polynomial isIts roots are $1$ and $1/7$, with the unit-modulus root simple, so it satisfies the root condition for a multistep method and is zero-stable. Combining this with the defect estimate proves convergence of order three, assuming sufficiently smooth $f$, order-three starting values, and the nearby branch of each implicit update.

More explicitly, the right-hand side is $h{F}_h(y_{n+2})$, where $F_h=6f/7-2h{f}^{\prime }f/7$ has a uniform local Lipschitz continuity bound for small $h$. A bounded zero-stable impulse response and the discrete Gronwall inequality therefore control the accumulated $O(h^4)$ defects by $O(h^3)$ over fixed time intervals. This is convergence of a zero-stable multiderivative method, rather than a direct application of the first-derivative-only Dahlquist equivalence theorem.

For the Dahlquist test equation, put $z=h\lambda$. The amplification roots satisfyWe show that no root can reach the unit circle for $\mathrm{Re}z<0$. If $\xi =e^{i\theta }$, solving the quadratic in $z$ givesWrite $c=\cos \theta$, $u=\mathrm{Re}C=-3/4+4c-c^2$, and $v=\mathrm{Im}C=(c-4)\sin \theta$. Then $u\le 9/4$ andSince $9/2-u>0$, this gives $|C|+u\le 9/2$. Consequently $|\mathrm{Re}\sqrt{C}{|}^2=(|C|+u)/2\le 9/4$, so both possible $z$ have nonnegative real part. Equality can occur only at $\theta =0$, which gives $z=0$ or $z=3$.

The leading coefficient has zeros $(3\pm i\sqrt{5})/2$, both in the right half-plane. Hence the amplification roots vary continuously as a pair throughout the left half-plane. At $z=-1$ they are $1/3$ and $1/5$; neither can leave the unit disk without crossing the unit circle, which the preceding calculation excludes. The same conclusion extends to the imaginary axis, with strict inequality away from $z=0$. At zero the roots $1,1/7$ satisfy the simplicity requirement.

Thus the method is A-stable. Its third order does not contradict the usual second-order barrier, because it is a multiderivative multistep method.

The tableau is the three-stage Lobatto IIIA method. For its stage matrix $A$, weights $b$, nodes $c$ and $C=\mathrm{diag}(c)$, the Butcher order conditions through order four arewhere $c^{\circ j}$ means coordinatewise powers. Direct substitution verifies all eight identities. Equivalently, this is a collocation Runge-Kutta method at the three Lobatto quadrature nodes, whose weights integrate cubic polynomials exactly.

Order five would require $b^T{c}^{\circ 4}=1/5$, but here $b^T{c}^{\circ 4}=5/24$. Therefore the method has exactly order four.

The stability function of a Runge-Kutta method is $R(z)=1+z{b}^T(I-zA{)}^{-1}1$. Substitution givesWriting the numerator and denominator as $N(z)$ and $D(z)$,The denominator has zeros $3\pm i\sqrt{3}$, both in the right half-plane. Therefore $|R(z)|\le 1$ throughout the closed left half-plane, with strict inequality in its interior. The method is A-stable. Since $R(z)\to 1$ at infinity, it is not L-stable.

For algebraic stability, the weights must be nonnegative and the real symmetric matrix $M=BA+A^{T}B-b{b}^T$, where $B=\mathrm{diag}(b)$, must be a positive semidefinite matrix.

The weights are positive, but $a_{11}=0$ and $b_1=1/6$, soA positive semidefinite matrix cannot have a negative diagonal entry, as its quadratic form at the first standard basis vector would be negative. Thus this Lobatto IIIA method is not algebraically stable, despite its A-stability.

Write $h=\Delta x$ and set $u_0=u_{M+1}=0$. In the discrete L2 norm $\Vert u{\Vert }_h^2=h{\sum }_{m=1}^M|u_m{|}^2$, the energy method givesTo obtain this identity, shift the indices in the centered second-difference sum. The centered first-difference matrix is skew-symmetric, so its contribution has zero real part even for complex grid values. The homogeneous Dirichlet boundary conditions remove the endpoint terms.

Therefore the semidiscretization is stable for every real $\alpha$ and every positive $h$, with the mesh-independent bound $\Vert u(t){\Vert }_h\le \Vert u(0){\Vert }_h$. No restriction on $|\alpha |h$ is needed for this energy estimate.

With $h=\Delta x$, $\tau =\Delta t$, $r=\tau /h^2$ and $c=\alpha \tau /h$, the Forward Euler method givesFor a Fourier mode $e^{im\theta }$, the amplification factor is $G(\theta )=1-4r{\sin }^2(\theta /2)+ic\sin \theta$. Set $s={\sin }^2(\theta /2)$. ThenThe expression in braces is an affine function of $s$, so it is nonpositive for all $s\in [0,1]$ precisely when its two endpoint values are nonpositive. For $\tau >0$, these conditions are $c^2\le 2r$ and $r\le 1/2$. Thus Forward Euler stability for centered advection-diffusion requires exactlyFor $\alpha =0$, only the first condition remains. If either bound fails, a band of Fourier modes has $|G|>1$, giving instability under von Neumann stability analysis. At equality the amplification bound still holds, so the endpoints are included. This is discrete L2 norm stability; nonnegative stencil weights would impose the additional spatial restriction $|\alpha |\Delta x\le 2$.

Keep the Courant number $\mu =\Delta t/\Delta x$ fixed and write $h=\Delta x$. An exact solution of $u_t=u_x$ is $u(x,t)=g(x+t)$. Substituting it into the update gives the defectwhere $s=x+t$ and $g^{(3)}$ is evaluated at $s$. The constant, first-derivative and second-derivative terms cancel. Dividing by $\Delta t=\mu h$ gives a second-order consistency error. Hence for fixed $0<\mu <1$, $\mu \ne 1/2$, the scheme has order two.

At the special value $\mu =1/2$, the update reduces to $u_m^{n+1}=u_{m+1}^{n-1}$. This is exactly the transport over two time steps, since $2\Delta t=\Delta x$, so the defect vanishes identically. With exact data at both starting levels it transports their values exactly; errors in those starting levels are still present, and the even and odd time subsequences evolve separately.

For a Fourier mode, the amplification polynomial isPut $a=1-2\mu$ and $\xi =e^{i\theta /2}\eta$. The roots becomeIf $|a|<1$, the square root is real, both roots have modulus one, and their separation is at least $2\sqrt{1-a^2}>0$. The companion matrix has eigenvectors $({\xi }_{\pm },1{)}^T$, so their uniform separation gives a mesh-independent bound on its powers. By power boundedness of a two-level Fourier scheme and the discrete Fourier transform, this proves stability for $0<\mu <1$.

If $|a|>1$, choose $\theta =\pi$. The roots have product of modulus one and unequal moduli, so one has modulus greater than one. The scheme is unstable.

The boundary cases require more than a root-modulus check. At $\theta =\pi$, $\mu =0$ gives $(\xi -1{)}^2=0$, while $\mu =1$ gives $(\xi +1{)}^2=0$. The companion matrices have nontrivial Jordan blocks and solutions of the form $(A+Bn)(\pm 1{)}^n$. Their linear growth rules out a uniform stability bound. This mode occurs on even periodic grids, which are enough to disprove mesh-uniform stability at the endpoints.

Therefore the stable range is exactlyWithin this range, second-order convergence additionally requires two appropriately accurate starting time levels. The unit-modulus parasitic root near $\theta =0$ does not decay, so poor initialization cannot be repaired by the scheme.

Let $x_m=mh$ be a uniform periodic mesh, with indices interpreted modulo the number of nodes. Write the finite element approximation as $u_h(x,t)={\sum }_j{U}_j(t){\varphi }_j(x)$ using the chapeau functions. The Galerkin method requiresThe mass matrix has $M_{mm}=2h/3$ and $M_{m,m\pm 1}=h/6$. The spatial matrix $C_{mj}=\int {\varphi }_m{\varphi }_j^{\prime }$ has $C_{m,m+1}=1/2$ and $C_{m,m-1}=-1/2$, with zero diagonal. Therefore the semidiscrete equations areor equivalently $M\dot{U}=CU$. The coefficients come from integrating the products of the two overlapping piecewise-linear hat functions on each interval; periodicity supplies the wraparound entries.

The mass matrix is a positive-definite matrix, and integration by parts with periodic boundary conditions gives $C^T=-C$. HenceThis is energy conservation for semidiscrete Galerkin advection: $U^{T}MU=\Vert u_h{\Vert }_{L^2(0,1)}^2$. Thus the semidiscretization is stable and conserves its finite element L2 norm.

For a mesh-independent comparison with the grid norm, the Fourier symbol of $M$ is $h(2+\cos \theta )/3$, between $h/3$ and $h$. ConsequentlyThe conserved energy therefore gives a uniform bound in $\sqrt{h{\sum }_m|U_m{|}^2}$ as well. Equivalently, each Fourier mode evolves with the purely imaginary exponent $\sigma (\theta )=3i\sin \theta /[h(2+\cos \theta )]$. This concerns the semidiscrete method; a chosen time integrator must separately control those imaginary eigenvalues.

A Runge-Kutta method with stage matrix $A$, weights $b$ and step size $h$ isIts order measures accuracy as $h\to 0$, while stability concerns how errors or stiff components propagate. Linear A-stability, strong stiff damping through L-stability, and nonlinear B-stability address different questions.

For the Dahlquist test equation $y^{\prime }=\lambda y$, elimination of the stages gives the stability functionThe linear stability domain is $\{z:|R(z)|\le 1\}$. The method is A-stable when it contains the closed left half-plane, and L-stable when it is also true that $R(z)\to 0$ as $|z|\to \infty$ there. These conditions describe, respectively, stability for every non-growing scalar linear mode and damping of very rapidly decaying modes.

The Forward Euler method has $R(z)=1+z$, so it fails A-stability, for example at $z=-3$. In fact no nontrivial explicit Runge-Kutta method is A-stable: its stability function is a nonconstant polynomial, which is unbounded along the negative real axis. The Backward Euler method has $R(z)=1/(1-z)$. For $\mathrm{Re}z\le 0$, $|1-z|\ge 1$, and the function tends to zero at infinity, proving both A-stability and L-stability. The implicit midpoint rule has $R(z)=(1+z/2)/(1-z/2)$. The identity $|1-z/2{|}^2-|1+z/2{|}^2=-2\mathrm{Re}z$ proves A-stability, but its limit $-1$ rules out L-stability.

For nonlinear equations, a dissipative vector field satisfies $⟨f(x)-f(y),x-y⟩\le 0$. The corresponding exact solutions contract because the time derivative of their squared distance is twice that inner product. A B-stable method preserves this contraction for every $h>0$, whenever its stages are well defined.

Let $d=y_n-{\tilde{y}}_n$, $D_i=Y_i-{\tilde{Y}}_i$, and $F_i=f(Y_i)-f({\tilde{Y}}_i)$. Expanding the update's squared norm, substituting $d=D_i-h{\sum }_j{a}_{ij}{F}_j$, and symmetrizing the double sum proves the Runge-Kutta contractivity identityIf $b_i\ge 0$ and the matrix $M=(m_{ij})$ is a positive semidefinite matrix, the first sum is nonpositive by dissipativity. The second double sum is nonnegative: sum the quadratic form of $M$ over each coordinate of the $F_i$. Thus algebraic stability implies B-stability. Applied to complex scalar linear equations with the real part of the Hermitian inner product, the same argument also implies A-stability, where the stages are defined.

The Backward Euler method has $b_1=a_{11}=1$, hence $M=(1)$, and is algebraically stable. One can prove its nonlinear contraction directly: $d_{n+1}=d+h(f(y_{n+1})-f({\tilde{y}}_{n+1}))$ gives $\Vert d_{n+1}{\Vert }^2\le ⟨d,d_{n+1}⟩\le \Vert d\Vert \Vert d_{n+1}\Vert$. The implicit midpoint rule has $a_{11}=1/2$, $b_1=1$ and $M=0$, so it is B-stable even though it is not L-stable.

Linear A-stability does not by itself imply nonlinear B-stability. The ODE trapezoidal rule has the same stability function as the implicit midpoint rule, but take the dissipative scalar field $f(y)=-\max (y,0)$ and $h=6$. For any positive starting value, its trapezoidal update is $y_{n+1}=y_n+3(f(y_n)+f(y_{n+1}))=-2y_n<0$. Thus two positive starting values have their distance doubled, violating B-stability. Its algebraic stability matrix also has a negative first diagonal entry. The field is globally Lipschitz continuous, so this is a well-defined nonlinear counterexample.

These examples distinguish accuracy, scalar stiff-mode stability, nonlinear contraction, and rapid stiff-mode damping. For implicit methods, solvability of the stage equations remains an additional requirement; a formal stability function alone does not establish it for arbitrary nonlinear problems.

Consider a bounded Lipschitz domain $\Omega$ and homogeneous Dirichlet boundary conditions forAssume $A$ is a bounded real symmetric matrix field with ${\xi }^{T}A(x)\xi \ge a_0|\xi {|}^2$ uniformly, $c\ge 0$ is bounded, and $f$ defines a bounded linear functional on $V=H_0^1(\Omega )$, the zero-boundary Sobolev space. The weak formulation isThe Cauchy-Schwarz inequality gives continuity $|a(u,v)|\le M\Vert u{\Vert }_V\Vert v{\Vert }_V$. Uniform ellipticity and the Poincaré inequality give a coercive bilinear form, $a(v,v)\ge \alpha \Vert v{\Vert }_V^2$. The Lax-Milgram theorem therefore supplies a unique weak solution, with $\Vert u{\Vert }_V\le \Vert \ell {\Vert }_{V^{\prime }}/\alpha$. Its proof represents $a$ by a bounded linear operator $T$ using the Riesz representation theorem; coercivity gives an injective operator with closed image, and a zero orthogonal complement makes that image all of $V$. This yields existence, uniqueness and the displayed estimate.

For a conforming finite element space $V_h\subseteq V$, the Galerkin method finds $u_h\in V_h$ such that $a(u_h,v_h)=\ell (v_h)$ for every $v_h\in V_h$. With a basis ${\varphi }_1,\dots ,{\varphi }_N$, this is the linear system $KU=F$, where the stiffness matrix is $K_{ij}=a({\varphi }_j,{\varphi }_i)$ and $F_i=\ell ({\varphi }_i)$. Its positive-definite matrix property proves unique discrete solvability. Local support gives a sparse matrix through elementwise assembly.

For the symmetric problem, the Ritz method minimizes $J(v)=a(v,v)/2-\ell (v)$ over $V_h$. Differentiating in every direction $v_h$ gives precisely the Galerkin method equations. Conversely, if $u_h$ solves them, $J(u_h+w_h)-J(u_h)=a(w_h,w_h)/2$, proving it is the unique minimum. The same argument over $V$ gives $J(v)-J(u)=\Vert v-u{\Vert }_a^2/2$, where $\Vert v{\Vert }_a=\sqrt{a(v,v)}$ is the energy norm. Thus the Ritz method chooses the best trial function in that norm.

Subtracting the continuous and discrete equations gives Galerkin orthogonality, $a(u-u_h,v_h)=0$. For arbitrary $v_h\in V_h$ and $e=u-u_h$, it follows that $a(e,e)=a(e,u-v_h)$. Coercivity and continuity prove the Céa lemma:In the symmetric energy norm, the Cauchy-Schwarz inequality improves the constant to one. Equivalently, Galerkin orthogonality gives the Pythagorean identity $\Vert u-v_h{\Vert }_a^2=\Vert u-u_h{\Vert }_a^2+\Vert u_h-v_h{\Vert }_a^2$. These estimates establish stability and show that approximation properties of the trial space determine convergence.

For continuous piecewise-linear elements on a shape-regular mesh in dimensions at most three, a finite element interpolation estimate gives $\Vert u-I_{h}u{\Vert }_{H^1}\le Ch\Vert u{\Vert }_{H^2}$ when $u\in H^2(\Omega )$. Combining it with the Céa lemma yields $O(h)$ energy norm error. More generally, degree-$p$ elements give $O(h^p)$ error in $H^1$ under the corresponding regularity and approximation assumptions. Mere mesh refinement cannot supply an order whose required solution regularity is absent.

An additional order in the L2 norm follows under a dual elliptic regularity assumption. Solve $a(v,z)=(e,v{)}_{L^2}$ and assume $\Vert z{\Vert }_{H^2}\le C\Vert e{\Vert }_{L^2}$. The Aubin–Nitsche duality argument uses Galerkin orthogonality to obtainTherefore $\Vert e{\Vert }_{L^2}\le Ch\Vert e{\Vert }_{H^1}$, giving $O(h^2)$ for the piecewise-linear case with $u\in H^2$. This improvement explicitly requires regularity of the dual problem.

As a concrete example, solve $-u^{\prime \prime }=1$ on $(0,1)$ with zero endpoint values. On a uniform mesh of spacing $h$, the piecewise-linear hat functions give $K_{jj}=2/h$, $K_{j,j\pm 1}=-1/h$ and $F_j=h$. The discrete equations are $2U_j-U_{j-1}-U_{j+1}=h^2$, whose solution is $U_j=x_j(1-x_j)/2$. The exact solution is $u(x)=x(1-x)/2$, so $u_h$ is its piecewise-linear nodal interpolant. On each interval $[x_j,x_{j+1}]$,Integrating the squared error and squared derivative error gives the explicit rateswhere the $H_0^1$ norm here is $\Vert v^{\prime }{\Vert }_{L^2}$.

For nonsymmetric coercive problems, the Galerkin method, Galerkin orthogonality and the Céa lemma still apply, but the symmetric quadratic minimization interpretation of the Ritz method generally does not. Nonconforming trial spaces, inexact integration and noncoercive equations need additional arguments beyond the conforming coercive theory proved here.