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

Insert the exact solution and expand about $t=t_{n+2}$. Since $f(y)=y^{\prime }$ and $f^{\prime }(y)f(y)=y^{\prime \prime }$, the left side isThe coefficient of $h^j{y}^{(j)}$ on the right side isfor $j\ge 1$, with the last term absent for $j=0$. These coefficients agree through $j=3$; at $j=4$ the left side minus the right side is $1/10$. Thus the local truncation error is $h^4{y}^{(4)}/10+O(h^5)$ and the method has order three.

At $h=0$ the first characteristic polynomial isIts roots satisfy the root condition for a multistep method: the only unit-modulus root is the simple root $1$, and the other root has modulus $1/5$. The method is consistent because part a gives positive order. The Dahlquist equivalence theorem, in the form allowed by the question for this derivative-augmented method, therefore proves convergence.

For the test equation $y^{\prime }=\lambda y$, put $z=h\lambda$. The amplification roots satisfyThe quadratic Schur stability criterion shows that both roots lie in the closed unit disk whenTo check these inequalities on the closed left half-plane, write $z=-u+iv$ with $u\ge 0$. The right side of the second inequality isAfter squaring, the difference between its square and $16|5\overline{D}+1+z{|}^2$ is a polynomial in $u$ and $v^2$ with nonnegative coefficients:The first inequality follows from the same displayed positive expression. The unit-circle cases are semisimple, so the entire closed left half-plane is in the linear stability domain. Hence the method is A-stable.

The Lagrange interpolation polynomial basis at $c_1,c_2$ isThe Collocation Runge-Kutta method has stages and updatewhere $a_{ij}={\int }_0^{c_i}{\ell }_j(s)ds$ and $b_j={\int }_0^1{\ell }_j(s)ds$. Explicitly,

A two-node collocation method already has order at least two. It has order at least three exactly when its node polynomial is orthogonal to constants on $[0,1]$:Equivalently,This is the additional quadrature order condition ${\sum }_i{b}_i{c}_i^2=1/3$.

Put $c_1=1/2$ and $c_2=c\ne 1/2$. The weights and stage matrix becomeFor algebraic stability of a Runge-Kutta method,soA positive-semidefinite matrix with a zero diagonal entry must have every entry in that row and column equal to zero. Here the off-diagonal entry never vanishes for finite $c\ne 1/2$. Therefore there is no admissible value of $c_2$ for which the method is algebraically stable.

This is the Backward Euler diffusion scheme with a centered second spatial difference. Taylor expansion gives first order in time and second order in space:With the parabolic scaling $\Delta t=\mu \Delta x^2$ and fixed $\mu$, the combined error is $O(\Delta x^2)$.

Insert the Fourier mode $u_m^n=G^n{e}^{im\theta }$ into the scheme. Von Neumann stability analysis givesFor every physical Courant number $\mu \ge 0$, one has $0<G(\theta )\le 1$. Thus the method is unconditionally stable: the range is $\boxed{\mu \ge 0}$.

Let there be $J$ interior points. The Dirichlet discrete Laplacian has positive valuesfor the eigenvalues of its negative, and the amplification eigenvalues are $(1+\mu {\eta }_j{)}^{-1}$. Therefore the Backward Euler diffusion stability on a finite Dirichlet interval isIf, as usual, a Courant number is restricted to nonnegative values, this again reduces to all $\mu \ge 0$. The additional negative branch is a finite-grid artefact and disappears to $-\infty$ as the mesh is refined.

The clamped energy space isTwice applying integration by parts, with the boundary terms killed by the clamped conditions, gives the symmetric bilinear formConsequentlyfor every nonzero $u\in \mathcal{H}$: equality forces $u^{\prime \prime }=0$ almost everywhere, and the clamped boundary values then force $u=0$. Hence $L$ is symmetric and positive definite.

Strictly under the stated assumption $p\in L^2$ rather than $L^{\infty }$, the first integral need not be finite for every $u\in H_0^2$. The literal energy domain is therefore $\{u\in H_0^2:\sqrt{p}{u}^{\prime \prime }\in L^2\}$; under the usual coefficient assumption $0<p_0\le p\le p_1<\infty$, it is exactly $H_0^2$ and the form is coercive there.

Define the energyIts first variation in direction $w\in \mathcal{H}$ isThus its minimizer $u$ satisfies the weak formulationwhich is the weak equation $Lu=f$. Positive definiteness makes $J$ strictly convex, so this stationary point is the unique minimizer; under uniform positivity of $p$, existence follows from the Lax-Milgram theorem.

Choose a mesh of $[-1,1]$ and the conforming space of Cubic Hermite finite elements: piecewise cubic functions that are globally $C^1$ and satisfy $v(\pm 1)=v^{\prime }(\pm 1)=0$. Let ${\varphi }_1,\dots ,{\varphi }_n$ be its nodal value-and-slope basis. For $u_n={\sum }_k{a}_k{\varphi }_k$, the Ritz method imposesHence the coefficient vector solvesThe stiffness matrix is symmetric positive definite by part a.

Write $H=-{\partial }_x^2+V(x)$. The real potential and periodic boundary conditions make $H$ a self-adjoint operator. Since $u_t=-iHu$,because the expectation of a self-adjoint operator is real. Equivalently, integrating the kinetic term by parts leaves $\int |u_x{|}^{2}dx$ and the periodic boundary term cancels. Thus the continuous Schrodinger equation preserves its $L^2$ norm.

The centered spatial second difference has error $O(\Delta x^2)$. The symmetric compositionis Strang splitting, so its local splitting error is $O(\Delta t^3)$ and its global time order is two. Because both semidiscrete subproblems are solved exactly, the total global error is

The periodic centered-difference kinetic matrix $T_h$ is real symmetric, and the sampled real potential $V_h$ is a real diagonal matrix. Therefore $-i{T}_h$ and $-i{V}_h$ are skew-Hermitian, so each matrix exponential in the Strang splitting is a unitary matrix. Their product is unitary as well. Hence

After spatial discretization, a linear time-dependent PDE produces an updateStability over $0\le n\Delta t\le T$ requires a bound $\Vert Q_h^n\Vert \le C_T$ independent of the mesh. If $Q_h$ is normal, the spectral theorem for normal operators givesso eigenvalue analysis is decisive. More generally, if $Q_h=X_h{\Lambda }_h{X}_h^{-1}$, thenThus eigenvalues suffice only when the eigenvector condition numbers are uniformly bounded and unit-circle eigenvalues are semisimple. Defective or increasingly nonnormal matrices can have large powers even though every eigenvalue lies in the unit disk.

For a constant-coefficient scheme on the whole line or a periodic grid, the discrete Fourier transform diagonalizes translation-invariant difference operators. This is Von Neumann stability analysis: insert $u_m^n=G(\theta {)}^n{e}^{im\theta }$ and require $|G(\theta )|\le 1$. Its advantages are simplicity, sharp mesh restrictions, and direct identification of unstable wavelengths. Its limitations are boundaries, variable coefficients, nonlinearities, and nonnormality; frozen-coefficient Fourier analysis then gives at most local evidence.

As a successful implicit example, the Backward Euler diffusion scheme haswhere the periodic or homogeneous-Dirichlet discrete Laplacian $D_h$ is symmetric negative semidefinite. Its eigenvalues are $(1-\Delta t{\lambda }_j(D_h){)}^{-1}\in (0,1]$, proving unconditional discrete-$2$-norm stability.

For a failure, consider explicit upwinding for $u_t+a{u}_x=0$ on a finite inflow grid:Its lower-triangular nonnormal upwind amplification matrix has only the eigenvalue $1-\mu$, so eigenvalues alone incorrectly suggest stability for $0\le \mu \le 2$. For $1<\mu <2$, however, interior alternating data are amplified by $|1-2\mu |>1$ before the boundary is felt, and the matrix powers have no mesh-uniform bound. The correct range is $0\le \mu \le 1$. This example isolates the missing hypothesis: spectral radius does not control powers of a nonnormal matrix family.

Suppose a semidiscrete linear PDE is $u_t=(A+B)u$. Direct evaluation of $e^{h(A+B)}$ may be expensive, while $e^{hA}$ and $e^{hB}$ can be cheap, parallelizable, or exactly structure-preserving. Splitting is especially effective when $A$ and $B$ represent different spatial directions, kinetic and potential energy, diffusion and reaction, or linear and nonlinear pieces.

The Lie-Trotter splitting commutator error follows from multiplying the exponential series:It therefore has local error $O(h^2)$ and global order one. Reversing the factors changes the sign of the leading commutator. The symmetric average in composition form is Strang splitting,whose symmetry $S_{-h}=S_h^{-1}$ removes even powers from the modified generator. Its local error is $O(h^3)$, involving nested commutators such as $[A,[A,B]]$ and $[B,[A,B]]$, and its global order is two.

For multidimensional diffusion, directional splitting replaces one large elliptic solve by successive one-dimensional tridiagonal solves; alternating-direction implicit methods are standard examples. For the Schrodinger equation, take $A=i{\partial }_x^2$ and $B=-iV(x)$. The kinetic subflow is diagonal in Fourier space, the potential subflow is pointwise multiplication, and Strang splitting is both second order and exactly norm-preserving.

A symmetric method has even global order. If its leading modified-flow defect is $h^{2p+1}E$, a symmetric composition $S_{a_{1}h}\cdots S_{a_{s}h}$ raises the order when $\sum a_i=1$ and $\sum a_i^{2p+1}=0$, together with the required higher commutator conditions. In particular, the higher-order composition of a symmetric splittingis fourth order because $2\gamma +\delta =1$ and $2{\gamma }^3+{\delta }^3=0$. Its negative middle step is harmless for reversible unitary problems but can be ill-posed or strongly unstable for parabolic semigroups. Higher-order parabolic splittings therefore use more specialized complex coefficients, commutator corrections, or extrapolation.