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

Substitute the exact solution and expand every value about $t_n$. The coefficient of $h^q{y}^{(q)}(t_n)$ in the defect is zero for $q=0,\dots ,p$ exactly whenThese are the order conditions for a linear multistep method, specialized to the two nonzero derivative coefficients. The first failed identity determines the leading local truncation error.

For $s=1$, imposing order two givesThe method is therefore the trapezoidal ruleIts first characteristic polynomial is $\rho (\zeta )=\zeta -1$, so it satisfies the root condition for a multistep method.

For $s=2$, the four conditions through order three giveand henceHereso the root condition for a multistep method again holds. Both methods are consistent and zero-stable, and the Dahlquist equivalence theorem therefore proves that both are convergent.

A highest-order method in this family has order $s+1$. The Second Dahlquist barrier says that an irreducible A-stable multistep method has order at most two, so $s+1\le 2$ and hence $s=1$. Part b then leaves only the trapezoidal rule. Its amplification factor isand $|R(z)|\le 1$ whenever $\mathrm{Re}z\le 0$. Thus the trapezoidal rule is the unique highest-order A-stable method of the stated form, apart from representations containing removable common factors.

For the nodes $c_1=\alpha$ and $c_2=1$, the Lagrange basis isDirect integration givesand thereforeThe Collocation Runge-Kutta method also requiresThis exposes a sign error in the printed tableau: its lower-right entry is shown as $(-1+2\alpha )/[2(1-\alpha )]$. With $1-2\alpha$ in that position, the tableau is exactly the claimed collocation method. Taken literally, the printed weights satisfy $b^{T}1=\alpha /(1-\alpha )$, so the method is not even consistent unless $\alpha =1/2$ and cannot be a collocation method for general $\alpha$.

For the intended collocation weights, interpolation makes the associated quadrature rule exact for every polynomial of degree at most one, so the method has order at least two. It has order at least three precisely when the node polynomial is orthogonal to constants:ThusAt $\alpha =1/3$ this is the two-stage Radau IIA method. Since one node is fixed at the endpoint, no value of $\alpha$ makes the two-node quadrature exact through degree three, so order four cannot occur.

For completeness, the tableau exactly as printed has order zero when $\alpha \ne 1/2$, since $b^{T}1\ne 1$. At $\alpha =1/2$ the two signs coincide because the second weight vanishes, and the resulting method has order two.

For algebraic stability of a Runge-Kutta method, the weights must be nonnegative andmust be positive semidefinite. With the intended collocation weights,Positive semidefiniteness is therefore possible only at $\alpha =1/3$; substitution gives nonnegative weights and a positive-semidefinite $M$. Hence the intended family is algebraically stable exactly whenWith the sign printed in the paper, one instead obtainsEquality forces $\alpha =1/2$, where $M$ still has nonzero off-diagonal entries and is indefinite. The literal printed tableau is consequently algebraically stable for no value of $\alpha$.

Write $H={\partial }_x^2-V(x)$. The real potential and periodic boundary condition make $H$ a self-adjoint operator on the periodic domain. Since $u_t=-iHu$,because the inner product $⟨u,Hu⟩$ is real. Thus the Schrodinger equation preserves the $L^2$ norm.

Let $H_h$ be the periodic centered second-difference matrix minus the real diagonal matrix containing $V(x_m)$. It is a Hermitian matrix, so the semidiscrete system iswith a skew-Hermitian matrix generator. ConsequentlyIts exact propagator $e^{-it{H}_h}$ is a unitary matrix, and therefore the semidiscretization is stable in the discrete $2$-norm, uniformly for all times and mesh sizes.

Apply the implicit midpoint rule, equivalently the Crank--Nicolson method, to the semidiscrete equation:It has order two. Its amplification matrix is the Cayley transformBecause $H_h$ is Hermitian, $Q$ is unitary. Hence $\Vert {\mathbf{u}}^{n+1}{\Vert }_2=\Vert {\mathbf{u}}^n{\Vert }_2$ exactly.

A real linear operator $\mathcal{L}$ on an inner-product space is a positive-definite operator when it is self-adjoint andfor every nonzero $u$ in its domain. In the variational setting one normally requires the stronger uniform estimate $⟨\mathcal{L}u,u⟩\ge c\Vert u{\Vert }_V^2$ for some $c>0$, which is coercivity.

Choose the Sobolev space $V$ encoding the homogeneous essential boundary conditions, setafter the appropriate integration by parts, and define the energy functionalIts first variation is $J^{\prime }(u)v=a(u,v)-\ell (v)$, so its stationary points are exactly the solutions of the weak formulationIf $a$ is bounded, symmetric, and coercive and $\ell \in V^{\prime }$, the Lax-Milgram theorem supplies a unique weak solution $u$. Moreover, for every $w\in V$,unless $w=0$. Thus $J$ is strictly convex, and $u$ is its unique global minimizer. This proves existence and uniqueness of the minimizer and of the weak solution simultaneously.

Let $\Omega =(-1,1{)}^2$ and use the clamped energy space $V=H_0^2(\Omega )$, whose traces satisfy $u={\partial }_{n}u=0$ on $\partial \Omega$. Two applications of integration by parts giveIf equality holds, then $\Delta u=0$. The maximum principle for harmonic functions and the zero Dirichlet boundary condition imply $u=0$. The biharmonic operator is therefore positive definite. The same identity and conclusion hold for the simply supported conditions $u=\Delta u=0$; either standard interpretation of the paper's phrase “zero boundary conditions” gives the result.

Near $t=0$, take $\Omega (t)=\log F(t)$ to be the continuous local matrix logarithm with $\Omega (0)=0$. The symmetry identity gives $F(-t)=F(t{)}^{-1}$, and uniqueness of this logarithm yieldsThus $\Omega$ is an odd function. The local-logarithm qualification is necessary because the matrix exponential is not globally injective; the statement is naturally understood either near $t=0$ or as an identity of formal power series.

For the Strang splittingreversing $t$ reverses all three factors, soIt is therefore a time-symmetric numerical method. Multiplication of the three exponential series shows agreement with $e^{t(A+B)}$ through degree two. Equivalently, the Baker--Campbell--Hausdorff formula gives an odd modified generatorfor a matrix $D$ made from nested commutators. Exponentiating givesfor a matrix $C$ depending on $A$ and $B$.

The three substeps have total length $2\alpha +(1-2\alpha )=1$. Their leading cubic defects add, so cancellation requiresTaking the real cube root and solving givesThis is the coefficient in the higher-order composition of a symmetric splitting; the middle substep $1-2\alpha$ is negative.

The composition is palindromic, and henceThus $G_{\alpha }$ is symmetric. Part c cancels the cubic term in its odd formal logarithm; symmetry forbids a fourth-degree term, so the next possible defect has degree five. Thereforewhich makes the composition a fourth-order method.

An $s$-stage Runge-Kutta method has stages and updateOn the Dahlquist test equation $y^{\prime }=\lambda y$, elimination of the stages gives the stability functionThe linear stability domain is the set where $|R(z)|\le 1$. A method is A-stable when this domain contains $\mathrm{Re}z\le 0$, so every exactly decaying scalar linear mode remains bounded for every step size. It is L-stable when it is A-stable and $R(z)\to 0$ as $|z|\to \infty$ in the left half-plane; this extra limit strongly damps unresolved stiff modes.

The rational function $R$ makes several useful conclusions immediate. No explicit Runge--Kutta method is A-stable because its stability function is a nonconstant polynomial and is therefore unbounded on the negative real axis. The implicit midpoint rule has $R(z)=(1+z/2)/(1-z/2)$ and is A-stable, but $R(z)\to -1$, so it is not L-stable. The Backward Euler method has $R(z)=(1-z{)}^{-1}$ and is L-stable. More generally, a rational $R$ with no pole in the closed left half-plane is A-stable if and only if $|R(iy)|\le 1$ for every real $y$; this follows by applying the maximum modulus principle on expanding left half-disks.

Scalar linear stability does not by itself control nonlinear perturbations. Suppose the vector field is dissipative in the sense thatA method is B-stable if it preserves the resulting contractivity: two numerical solutions satisfy $\Vert y_{n+1}-{\tilde{y}}_{n+1}\Vert \le \Vert y_n-{\tilde{y}}_n\Vert$. A practical sufficient condition is algebraic stability of a Runge-Kutta method: $b_i\ge 0$ andTo prove the implication, let $D_i=Y_i-{\tilde{Y}}_i$ and $F_i=f(Y_i)-f({\tilde{Y}}_i)$. Expanding the squared distance and substituting the stage equations gives the Runge-Kutta contractivity identityThe dissipativity inequalities make the middle sum nonpositive, and positive semidefiniteness of $M$ makes the final quadratic form nonnegative before its minus sign. The distance therefore cannot increase. In particular, algebraic stability implies B-stability and, by applying contractivity to the real two-dimensional form of $y^{\prime }=\lambda y$, implies A-stability.

Important collocation families illustrate these notions. Gauss methods are A-stable, symmetric, and have order $2s$, but they do not damp infinitely stiff modes. Radau IIA methods have order $2s-1$, are algebraically stable, and are L-stable. These properties explain why A-stability controls unrestricted linear decay, L-stability is useful for stiff transients, and algebraic or B-stability is the stronger tool for nonlinear dissipative equations.

Consider first a scalar constant-coefficient Cauchy problem for a partial differential equation on the whole line. A translation-invariant spatial discretization is a convolution operator, so the discrete Fourier transform turns it into multiplication by its Fourier symbol. This reduction is exact because Fourier modes are the simultaneous generalized eigenfunctions of every translation-invariant stencil.

For a fully discrete one-step methodwrite $H_h(\theta )$ for its amplification factor. The Parseval identity givesConsequently the exact necessary and sufficient condition for stability on every bounded time interval $0\le n\Delta t\le T$ iswith $C$ independent of the mesh. Sufficiency follows directly from Parseval; necessity follows by choosing transformed initial data concentrated where the multiplier is largest. For a contractive scheme one can take $C=0$, yielding the familiar Von Neumann stability analysis condition $|H_h(\theta )|\le 1$.

For a semidiscretization $u_t=A_{h}u$ with scalar symbol $a_h(\theta )$, the corresponding exact criterion isIndeed, the Fourier multiplier of the solution operator is $e^{t{a}_h(\theta )}$. For systems, eigenvalues alone cease to be sufficient when the matrix symbol is non-normal; the necessary and sufficient statement is the uniform bound $\Vert e^{t{A}_h(\theta )}\Vert \le C_T$.

On the whole lattice, a finite stencil defines a Laurent operator, whose operator norm is the essential supremum of its symbol. Restricting the same stencil to a half-line gives a Toeplitz operator. The Cauchy symbol condition remains necessary for a stable initial-boundary scheme because data localized far from the boundary behave like the whole-line problem for a finite time. It is not sufficient: the boundary closure can support growing modes or amplify incoming modes even when every interior Fourier mode is stable.

Boundary stability therefore requires a uniform estimate for the forced half-line recurrence. After a Laplace transform in time and a Fourier transform in tangential variables, one solves a normal-direction recurrence. The Uniform Kreiss--Lopatinskii condition requires the boundary equations to determine its decaying roots with a uniformly bounded inverse. In practical terms, one imposes one independent boundary condition for each incoming characteristic or numerical mode and none for outgoing modes. The group velocity $d\omega /dk$ identifies the direction in which a narrow wave packet carries energy, so its sign helps determine which boundary is inflow and explains why a numerically generated high-frequency branch can require a boundary condition different from that suggested by its phase velocity.

For a two-step method, the Fourier substitution produces an amplification polynomial of a multilevel finite difference scheme. Every root must lie in the closed unit disk, and unit-modulus roots must be simple, uniformly in the wavenumber. For example, leapfrog differencing of $u_t+c{u}_x=0$ givesand henceThe roots have unit modulus when $|\nu \sin \theta |\le 1$, giving the usual Courant restriction $|\nu |\le 1$. At the endpoint, a wavenumber with $|\nu \sin \theta |=1$ produces a repeated unit root and violates the uniform root condition; strict $|\nu |<1$ avoids this marginal linear growth. The second root is the familiar oscillatory computational mode, illustrating why a multilevel scheme requires its full amplification polynomial rather than a single multiplier.