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

Along an exact autonomous-ODE solution, $f(y)=y^{\prime }$ andMove every term to the left and substitute the Taylor expansions about $t_n$. The coefficients of $h^j{y}^{(j)}(t_n)$ vanish for $0\le j\le 4$, while the first nonzero coefficient isThus the local defect is $O(h^5)$. At $h=0$ the first characteristic polynomial iswhich satisfies the root condition for a multistep method. The method is therefore zero-stable and has

Apply the method to the Dahlquist test equation $y^{\prime }=\lambda y$ and set $z=h\lambda$. Its amplification roots satisfyAt $z=0$ the roots are $1$ and $5/11$. A root can leave the unit disk only through $\xi =e^{i\phi }$. Substitution gives the boundary equationDirect separation into real and imaginary parts shows that both branches satisfy $\mathrm{Re}z\ge 0$; equality occurs on the branch through $z=0$. Hence no root crosses the unit circle in $\mathrm{Re}z<0$. Moreover both roots tend to zero as $z\to \infty$ in the left half-plane. Consequently the method is A-stable, and it also strongly damps the infinitely stiff limit.

Let $A=(a_{ij})$, $b=(b_i)$, $1=(1,1,1{)}^T$, and $c=A1$. From the displayed tableau,The quadrature moments satisfyDirect substitution into the remaining Butcher order conditions gives, for example,and all rooted-tree conditions of orders at most five are satisfied. The order-six moment already fails:The Runge-Kutta method therefore has

Algebraic stability requires $b_i\ge 0$ and positive semidefiniteness ofThe three weights are positive, but direct calculation givesA positive-semidefinite matrix cannot have a negative diagonal entry. Hence

The equation can be writtenFor the squared L2 norm, periodic integration by parts and the reality of $V$ giveThus the Schrodinger equation generates a unitary flow and $\Vert u(\cdot ,t){\Vert }_2$ is constant.

Here $\Delta x=1/M$ and periodic indexing is taken modulo $2M$. The second-order central difference matrix is the real symmetric circulant matrixSince $V$ is real diagonal, $H=A-V$ is Hermitian. Therefore $iH$ is skew-Hermitian, and

Split the generator as $iA-iV$. One Strang splitting step isThe two potential half-steps are componentwise multiplications. The circulant matrix $A$ is diagonalized by the discrete Fourier transform; its eigenvalues areThus the middle step consists of a Fast Fourier transform, multiplication of Fourier coefficient $j$ by $e^{ih{\lambda }_j}$, and an inverse transform. Its cost is $O(M\log M)$ and every factor is unitary, so the implementation preserves the discrete norm exactly up to roundoff.

Let $L_h=D_{xx}+\alpha D_x$, where the centered second-difference matrix $D_{xx}$ is real symmetric negative definite under homogeneous Dirichlet boundary conditions, and the centered first-difference matrix $D_x$ is real skew-symmetric. HenceThe semidiscrete energy method gives $\Vert \mathbf{u}(t){\Vert }_2\le \Vert \mathbf{u}(0){\Vert }_2$ for every real $\alpha$, so the scheme is stable.

For the Cauchy problem, insert the Fourier mode $u_m=e^{im\vartheta }$. The spatial symbol isIts real part is nonpositive for every $\vartheta$, so each mode has modulus $e^{t\mathrm{Re}{\lambda }_h}\le 1$. By the discrete Fourier transform and Parseval identity, the scheme is stable in the discrete $L^2$ norm.

Forward Euler method giveswhereIts amplification factor isWriting $X={\sin }^2(\vartheta /2)$, the condition $|G{|}^2\le 1$ for every $0\le X\le 1$ is equivalent toThereforewith the second bound omitted when $\alpha =0$.

Write the operator in Sturm-Liouville form:The natural solution space is the Sobolev space $H_0^1(0,1)$. Multiplication by a test function $v\in H_0^1(0,1)$ and integration by parts gives the symmetric bilinear formIt is bounded andby the Poincare inequality. Thus $L$ is positive definite and $a$ is coercive. The Lax-Milgram theorem gives a unique weak solution of the variational problem

Let $\{{\varphi }_j{\}}_{j=1}^{N-1}$ be the piecewise-linear hat functions on a mesh, and putThe Ritz method requireswhere, for $f\equiv 1$,Local support makes $K$ symmetric tridiagonal, and coercivity makes it positive definite.

Stability asks whether perturbations already present in data, arithmetic, or an earlier numerical step remain controlled under subsequent evolution. It complements consistency of a numerical method: consistency says that an exact smooth solution nearly satisfies one numerical step, while stability prevents the accumulated local defects from being amplified without bound. For a well-posed differential equation these two properties are what make convergence possible.

For a linear multistep method, zero-stability concerns the limit $h\to 0$. The roots of its first characteristic polynomial must lie in the closed unit disk, and every root on the unit circle must be simple. The Dahlquist equivalence theorem states that a consistent linear multistep method converges precisely when it is zero-stable. The parasitic root $5/11$ in Question 1 decays, whereas a repeated root at $1$ would turn small defects into secular growth.

Absolute stability instead fixes $z=h\lambda$ for the Dahlquist test equation. A one-step method has amplification factor $R(z)$ and absolute-stability regionA-stability means that $\mathcal{S}$ contains the entire closed left half-plane, matching every decaying scalar linear problem. L-stability additionally requires $R(z)\to 0$ as $z\to -\infty$, which suppresses unresolved fast transients in a stiff differential equation. Backward Euler is L-stable; the trapezoidal rule is A-stable but approaches $-1$ and can retain stiff oscillations; forward Euler is stable only in the disk $|1+z|\le 1$.

For nonlinear dissipative systems, scalar absolute stability can be insufficient. B-stability controls distances between numerical solutions of contractive differential equations. For Runge–Kutta methods, algebraic stability—nonnegative weights and positive semidefiniteness of $BA+A^{T}B-b{b}^T$—is a useful sufficient condition for B-stability. Question 2 shows how a single negative diagonal entry disproves it.

After spatial discretization of a partial differential equation, stability can be studied through the semidiscrete matrix. If its Hermitian part is nonpositive, the energy method proves contractivity without diagonalizing it. A skew-Hermitian generator instead conserves norm, as in the discrete Schrodinger problem of Question 3. Eigenvalues alone can be misleading for a non-normal matrix: transient growth may be large even when every eigenvalue lies in the left half-plane, so matrix norms, logarithmic norms, resolvent bounds, or a pseudospectrum may be needed.

For constant-coefficient Cauchy problems, Von Neumann stability analysis inserts Fourier modes and bounds their amplification factors. Question 4 gives a typical result: centered diffusion contributes a negative real symbol and centered advection an imaginary symbol. Semidiscrete evolution is stable, yet forward Euler imposes both the parabolic restriction $\Delta t=O(\Delta x^2)$ and an advection-diffusion restriction. This illustrates a Courant–Friedrichs–Lewy condition: a stable spatial approximation need not remain stable under an arbitrary time stepper.

For a well-posed linear initial-value problem and a consistent finite-difference approximation, the Lax equivalence theorem identifies stability with convergence. Its hypotheses matter: it does not by itself cover nonlinear equations, inconsistent boundary closures, changing norms, or non-smooth solutions. Boundaries also defeat a naive whole-line Fourier argument; an energy estimate, normal-mode boundary analysis, or a discrete semigroup bound must include the boundary treatment.

Practical stability analysis therefore starts from the structure of the differential equation. Conservation laws favor unitary or symplectic methods, diffusion favors A- or L-stable implicit methods, monotone transport may require strong-stability-preserving time stepping and upwind fluxes, and stiff splitting requires attention to both the factors and their commutators. Stability does not guarantee accuracy: a heavily damped method may be stable while erasing the solution, and a stable computation with $h\lambda$ near the edge of its stability region may have an unacceptable phase error.