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

Introduce one stage $k=f(t_n+h/2,Y)$ with $Y=y_n+hk/2$, and update $y_{n+1}=y_n+hk$. This is the implicit midpoint rule as a one-stage Runge-Kutta method, with Butcher tableauFor a smooth vector field, write $\Delta =y_{n+1}-y_n$. Expanding its stage equation giveswhich matches the exact Taylor expansion through degree two. The method therefore has order at least two. On the Dahlquist test equation $y^{\prime }=\lambda y$, its stability function isThe third coefficient differs from the exponential's $1/6$, so the order is exactly two.

There is no pole in $\mathrm{Re}z\le 0$, andThus $|R(z)|\le 1$ throughout the closed left half-plane:This conclusion is A-stability, not strong damping of arbitrarily stiff modes: $R(z)\to -1$ along the negative real axis.

Let $Q(y)=y^{T}Sy$. In the usual pointwise meaning of an ODE invariant, differentiation along an exact solution through an arbitrary state givesFor a nonautonomous equation this identity must hold at the stage time as well as at the initial time. It follows directly from invariance for solutions started at any time and state, or from the stated initial-time invariance when the exact flow onto the stage-time domain is invertible. This is the first-integral condition used below.

Put $Y=(y_n+y_{n+1})/2$ and $t_{\ast }=t_n+h/2$. Since $S$ is symmetric,Induction over all well-defined implicit steps givesNo positive-definiteness assumption on $S$ was used. Nor is uniqueness of the stage essential to this algebraic identity: any stage solution satisfying the method and the pointwise invariant condition preserves $Q$. This is the one-stage case of Runge-Kutta conservation of quadratic invariants.

Use the coefficient $2-\alpha$ in the numerator, as verified in the PDF; the converted TeX lost its $-\alpha$. Multiplying the method by $2+\alpha$ gives the characteristic polynomials of a linear multistep methodFor every admissible real $\alpha$, $\rho (1)=0$ and ${\rho }^{\prime }(1)=2\alpha =\sigma (1)$, giving consistency of a numerical method. Factorization givesThe second root is $r=(2-\alpha )/(2+\alpha )$. The root condition for a multistep method requires $|r|\le 1$ and forbids a repeated root at one. For real $\alpha \ne -2$,At $\alpha =0$ both roots equal one, violating zero-stability. For $\alpha >0$ the second root lies strictly inside the unit disk. The Dahlquist equivalence theorem therefore givesTo find the order rather than only consistency, expand the exact-solution residual through the generating function:Thus its unnormalized local truncation error begins with $(\alpha /3)h^3{y}^{(3)}$ when $\alpha >0$, and the global method order is exactly two for every convergent parameter. For comparison, $\alpha =0$ has formal order three but is not convergent because of its double unit root:These orders assume a smooth solution, the usual local Lipschitz condition and starting approximations of the required accuracy.

Only $\alpha >0$ needs consideration, by the preceding zero-stability analysis. On the Dahlquist test equation, all roots of the amplification polynomialmust satisfy the unit-disk root condition if the method is A-stable. Evaluate it at $\zeta =-1$:Take a real negative $z<-4/\alpha$. Then $P_z(-1)<0$, whereas the leading coefficient $2+\alpha -z$ is positive and $P_z(\zeta )\to +\infty$ as $\zeta \to -\infty$. The intermediate value theorem gives a real root strictly below $-1$. Its numerical mode grows in modulus, contradicting absolute stability at this negative test parameter. ConsequentlyThis proof detects an unstable root directly; the second-order value alone would not exclude A-stability. Equivalently the limiting derivative polynomial has an exterior root $-\alpha -\sqrt{1+{\alpha }^2}$, explaining the instability for large negative $z$.

There is a genuine indexing defect in the printed PDF: its new value is at $m+1$, while every old-time coefficient is centered at $m$. On a fixed grid the literal residual containsFor a general solution the leading term is $u_x/(\mu \delta )$, so this is not a consistent diffusion discretization. It also leaves one new interior value unspecified while imposing a generally incompatible update on a prescribed boundary value. The calculation that follows is for the intended same-node update $U_m^{n+1}=U_m^n+k{L}_{\delta }{U}_m^n$; that correction is necessary, not an unnoticed transcription change.

For the midpoint flux stencil for one-dimensional diffusion, defineAssume sufficient smoothness of $a$ and the solution for the displayed remainders. Multiplying the centered Taylor expansions of $a$ and $u$ givesThe first time difference is $(u(t+k)-u(t))/k=u_t+k{u}_{tt}/2+O(k^2)$. After using the diffusion equation, the normalized local truncation error isTherefore, for the corrected update,The first is the residual divided by the time step; the second is the unnormalized one-step defect. At fixed $\mu =k/{\delta }^2$ they are respectively $O({\delta }^2)$ and $O({\delta }^4)$. A rate based on derivatives is not justified for a merely bounded nonsmooth $a$.

Use the corrected same-node method identified above, with $\delta =1/M$, interior indices $1,\dots ,M-1$ and zero endpoint values. Its spatial matrix $L_{\delta }$ is real symmetric. Summation by parts gives, for a vector extended by zero at the endpoints,Hence $-L_{\delta }$ is positive definite, and ${\sum }_j(v_{j+1}-v_j{)}^2\le 4{\sum }_j{v}_j^2$ bounds its eigenvalues by $4a_{+}/{\delta }^2$. The forward Euler method has amplification matrix $I+k{L}_{\delta }$. Its eigenvalues lie in $[1-4\mu a_{+},1]$, so their moduli are at most one whenSymmetry makes this a contraction in the mesh-weighted discrete L2 norm, uniformly over all grids and all coefficient samples satisfying the bounds. There is also a maximum norm proof: the update weights are $\mu a_{m-1/2}$, $1-\mu (a_{m-1/2}+a_{m+1/2})$, and $\mu a_{m+1/2}$. They are nonnegative and sum to one, giving a contraction after the boundary values are included.

This mesh-independent bound is sharp. For $a\equiv a_{+}$, the spatial eigenvectors are $\sin (j\pi m/M)$ and their amplification factors areIf $\mu >1/(2a_{+})$, sufficiently fine grids have a highest-frequency factor below $-1$, giving exponential growth in the step number.

On one fixed grid, rather than uniformly under refinement, the sharp universal stability of a numerical method bound in the discrete L2 norm is slightly larger:Indeed the quadratic form above is bounded by $a_{+}$ times the constant-coefficient difference form, whose largest eigenvalue is $4{\cos }^2(\pi /(2M))/{\delta }^2$; equality is realized by the constant coefficient. Its limit as $M\to \infty$ is the boxed bound. No such fixed-grid matrix conclusion repairs the ill-defined shifted-boundary update in the literal printed formula.

Put $\nu =k/\delta$ and use the discrete Fourier transform on the infinite grid. A spatial Fourier mode $e^{i(j\theta +\ell \varphi )}$ has time coefficient satisfyingIts amplification polynomial of a multilevel finite difference scheme is ${\zeta }^2-2is\zeta -1$, with rootsFor $|s|<1$ they are distinct and have modulus one. Since $|s|\le 2\nu$, a fixed $0<\nu <1/2$ bounds their separation below by $2\sqrt{1-4{\nu }^2}$. Solving for the two mode amplitudes therefore bounds every power of the two-level companion matrix uniformly in frequency and step number. Parseval identity transfers this bound to the discrete L2 norm.

If $\nu >1/2$, an open set of frequencies around $(\pi /2,\pi /2)$ has $|s|>1$ and a root of modulus greater than one. Fourier data supported there grow exponentially. At $\nu =1/2$, the central frequency has the double root $i$, and the recurrence admits $v_n=(A+Bn)i^n$. Thus the root moduli being one do not suffice: a nontrivial Jordan block makes the two-level powers grow linearly. Although this exact frequency has measure zero on the infinite grid, the matrices depend continuously on frequency. For every step count, a neighborhood has powers nearly as large, so their essential supremum is still unbounded. Wave packets in those neighborhoods rule out a uniform stability of a numerical method estimate in the discrete L2 norm.

For arbitrary perturbations in both starting levels the sharp positive range is thereforeThe zero-step limiting recurrence is bounded but does not advance time. A specially filtered startup can suppress a repeated-root mode; that restriction does not establish stability of the full two-level scheme. This is the two-dimensional leapfrog stability threshold.

Divide the recurrence by $2k$ and substitute a smooth exact solution. Centered Taylor expansions in time and space giveThus the leapfrog advection scheme is consistent with normalized residual $O(k^2+{\delta }^2)$. Its unnormalized defect is $O(k^3+k{\delta }^2)$. For a fixed $0<\nu <1/2$, the uniform two-level bound proved above and accumulation over $O(1/k)$ steps yield, for a sufficiently smooth solution,Both starting levels must be supplied consistently; errors of order $k^2+{\delta }^2$ give second-order convergence. For general square-integrable functions as initial data, point sampling is not defined on equivalence classes. Use cell averages or another bounded grid projection, together with a bounded consistent startup. Approximate the initial data by smooth functions, apply the smooth-data convergence estimate, and use the uniform numerical and exact translation bounds to pass to the L2 space limit. This is the content of Lax equivalence theorem in the present setting.

At $\nu >1/2$ the growing Fourier mode prevents general convergence. At the endpoint, initial data supported in increasingly narrow neighborhoods of the double root can have a second-level error of norm $O(k)$ but a final error of order one after $O(1/k)$ steps. These are vanishing starting errors, so the endpoint fails general two-level convergence as well. Consequentlyis also the general convergence range, with a consistent stable initialization. Claims for a particular specially chosen endpoint startup would be a different, restricted assertion.

Assume $n\ge 3$ and use periodic indices for the Jacobi method. At iteration $r$, it isAll values on its right side belong to the old iterate. For lexicographic Gauss-Seidel iteration, compute in the order $1,\dots ,n$ and use a newly available coordinate immediately:The wrap-around term in the last row uses the new $x_1$, while the wrap-around term in the first row uses the old $x_n$. This is important for the actual iteration matrix and its convergence behavior.

The circulant matrix is singular, so the meaning of convergence needs care. If $C$ is the cyclic shift, $A=C+C^{-1}-2I$ has Fourier eigenvalues $-4{\sin }^2(\pi j/n)$. Its null space consists of constants. Therefore $Ax=b$ is solvable exactly when ${\sum }_j{b}_j=0$, and any solution is determined only up to a constant. Neither iteration can converge to a unique solution for arbitrary right-hand sides.

For a compatible $b$, fix one solution $x_{\ast }$ and examine the error. The iteration matrix of the Jacobi method is $(C+C^{-1})/2$, with eigenvalues $\cos (2\pi j/n)$. The constant mode has eigenvalue one. For odd $n$ all other eigenvalues have modulus less than one, so the iterates converge to $x_{\ast }+c1$, with the constant selected by the initial mean. For even $n$, the alternating vector has eigenvalue $-1$, producing an undamped oscillation. With $b=0$ and alternating initial vector this is an explicit nonconvergent example. For even $n$ convergence occurs only when that error component vanishes. When $\sum b_j\ne 0$, the mean instead drifts by $-(\sum b_j)/(2n)$ at each step, so Jacobi cannot converge.

For Gauss-Seidel iteration, put $K=-A=D+L+U$, so $K$ is positive semidefinite with diagonal two. The homogeneous error iteration is $e^{+}=Be$, $B=-(D+L{)}^{-1}U$. Updating a coordinate exactly minimizes $E(e)=e^{T}Ke/2$ in that coordinate; completing the square gives an energy decrement $|\Delta e_j{|}^2$. Summing within a sweep givesThe same identity holds for the Hermitian energy of complex vectors. If $Be=\lambda e$ and $e$ is nonconstant, its energy is positive, so $|\lambda |\le 1$. Equality forces every update to be zero, whence $Ke=0$ and $e$ is constant, a contradiction. Thus all nonconstant eigenmodes have modulus strictly below one.

The eigenvalue one has only the constant eigenspace. To exclude a generalized eigenvector at one, observe thatIf $(B-I)v=1$, multiplication by $w^T$ would give $0=n$. Hence the constant eigenvalue is algebraically simple, andFor compatible $b$, Gauss-Seidel therefore converges from every start to some solution, for either parity of $n$. Its limit differs from $x_{\ast }$ by $1w^T(x^{(0)}-x_{\ast })/n$. For an incompatible $b$, the affine update gives $w^T{x}^{(r+1)}-w^T{x}^{(r)}=-{\sum }_j{b}_j$, so the iterates again cannot converge.

Thus semiconvergence of cyclic Poisson iterations gives the precise answer:Both have spectral radius one on the unrestricted space, so neither meets the strict stationary-iteration criterion for a unique limit independent of the starting constant. A spectral radius of one alone is not a proof that Gauss-Seidel's compatible iterates fail to converge. A normalization or projection removing the constant mode makes this distinction explicit.

Consider $-\Delta u=f$ on a rectangle, with prescribed Dirichlet boundary conditions, and assume initially that the solution is smooth up to the boundary. A finite difference method chooses grid values as unknowns, approximates the differential operator by a local stencil, incorporates the boundary data, and solves the resulting sparse linear system. Accuracy of a stencil by itself is insufficient: one also needs a well-posed discrete problem, a stability estimate uniform under refinement, and a sufficiently accurate algebraic solve.

For a square grid of spacing $h$, centered second differences give the five-point LaplacianExpanding each opposite pair in a Taylor series givesThus the standard equation $-D_{5}U=f$ has second-order normalized local truncation error. Known boundary values move to the right-hand side; they are not iterated as independent unknowns. Nonrectangular geometry requires a consistent boundary approximation rather than silently retaining a full interior stencil outside the domain.

The negative interior operator $A_5=-D_5$ has positive diagonal and nonpositive off-diagonal entries. For a grid vector extended by zero on the boundary,The edge sum includes links to the fixed boundary. Connectivity makes the matrix positive definite, proving uniqueness and providing an energy structure. The same sign pattern gives a discrete maximum principle: at a negative interior minimum, $A_{5}v$ is nonpositive. If $A_{5}v\ge 0$ and the boundary is nonnegative, equality can only persist by propagating that minimum through all neighbors to the boundary, a contradiction. Therefore $v\ge 0$.

On the unit square, the barrier $b(x,y)=[x(1-x)+y(1-y)]/4$ satisfies $A_{5}b=1$ exactly. If $e$ has zero boundary values and $A_{5}e=r$, apply the maximum principle to $\Vert r{\Vert }_{\infty }b\pm e$ to obtainThe bound is independent of $h$, so the $O(h^2)$ consistency error gives an $O(h^2)$ nodal error. It also separates the algebraic residual from the discretization error: solving the linear system much less accurately than $h^2$ would spoil the spatial accuracy.

A direct higher-order design is possible. In one coordinate,has fourth-order accuracy. Adding the corresponding formula in the other coordinate gives a wider stencil. Its advantages must be balanced against boundary closures, more distant couplings, and the possible loss of monotonicity. Higher order does not automatically imply a stable or geometrically convenient scheme.

The Mehrstellen method takes a different route: it uses the Poisson equation to eliminate leading truncation derivatives while keeping a compact unknown stencil. Define the nine-point finite-difference stencilHere each sum has four terms. Equivalently $D_9=D_5+(h^2/6){\delta }_{xx}{\delta }_{yy}$. Expanding givesThe uncorrected nine-point operator is only second order for a general function. But $-\Delta u=f$ implies ${\Delta }^{2}u=-\Delta f$. Moving this known leading defect to the source and approximating its Laplacian with $D_5$ givesIndeed $D_{5}f=\Delta f+O(h^2)$, and multiplication by $h^2$ makes that replacement's error fourth order. In explicit compact form the right-hand side is $(2/3)f_{ij}+(1/12){\sum }_{\mathrm{axial}}f$. Only known source samples, not additional distant solution unknowns, enter this correction. The corresponding unscaled stencil residual is $O(h^6)$. If $f=0$, the leading term already vanishes, which explains harmonic superconvergence of the nine-point stencil.

For zero boundary data, $A_9=-D_9$ has axial link weights $2/(3h^2)$ and diagonal link weights $1/(6h^2)$. Their energy sum proves positive definiteness, and their sign pattern again proves the discrete maximum principle. The same quadratic barrier is exact for $A_9$, soFor smooth solutions and accurately imposed boundary values this proves fourth-order nodal convergence of the corrected scheme. Singular corners, rough sources or low-order boundary treatment can reduce the observed rate; a formal interior expansion is not a global regularity theorem.

More generally, compact high-order schemes choose a solution stencil and a source stencil together, match their Taylor expansions, and use differentiated forms of the PDE to remove derivatives that would otherwise need a wide solution stencil. Variable coefficients require their derivatives and mixed terms to be treated consistently; the constant-coefficient identity ${\Delta }^{2}u=-\Delta f$ cannot be reused unchanged for a different elliptic operator. On irregular meshes, conservation and consistent fluxes may be more useful starting points than a symmetric Cartesian formula.

Finally, discretization and solution cost are coupled. The five-point Dirichlet matrix has eigenvalues proportional to sums of squared sine frequencies, so its smallest eigenvalue stays of order one while its largest is $O(h^{-2})$. Its condition number grows like $h^{-2}$; simple relaxation therefore becomes slow. Sparse direct elimination, conjugate gradient method with appropriate matrix preconditioning, and the multigrid method are alternatives. A sound design combines local accuracy, boundary fidelity, mesh-independent stability, and an algebraic solver whose residual is small enough for the intended error tolerance.

The multigrid method addresses a characteristic failure of single-grid relaxation: its cheapest updates rapidly remove oscillatory error but barely change smooth error. For example, the one-dimensional Poisson equation matrix has diagonal $2/h^2$. The weighted Jacobi method with weight $\omega$ has Fourier error factorAt low physical frequency, $\theta =O(h)$ and this factor is $1-O(h^2)$, requiring $O(h^{-2})$ sweeps for a fixed error reduction. With $\omega =2/3$, high grid frequencies $\pi /2\le |\theta |\le \pi$ are multiplied by at most $1/3$ in modulus. Smooth error on a fine mesh is well represented on a coarser mesh, where it is relatively more oscillatory and hence cheaper to remove. These two complementary mechanisms motivate smoothing and coarse-grid correction.

Write the discrete problem as $A_h{u}_h=f_h$. For an approximation $v_h$, the residual is $r_h=f_h-A_h{v}_h$ and the error $e_h=u_h-v_h$ satisfies $A_h{e}_h=r_h$. It is this error equation that is transferred to a coarse grid; restricting the approximate solution alone would not solve the missing error. With prolongation $P$ and restriction $R$, solveFor nested Cartesian grids with $H=2h$, linear interpolation copies a coarse value at a coincident fine node and averages the two nearest coarse values at an intermediate node. Its two-dimensional counterpart is bilinear interpolation. Full-weighting restriction uses weights $(1,2,1)/4$ in one dimension andin two dimensions. With these choices $R=2^{-d}{P}^T$ in ordinary coordinate inner products. A Galerkin coarse-grid operator is $A_H=R{A}_{h}P$; alternatively one can rediscretize the PDE on the coarse mesh when that choice remains compatible with the transfers.

For a symmetric positive-definite fine operator and $R=c{P}^T$, the exact correction matrix isIt satisfies $C_{h}P=0$ and $P^T{A}_h{C}_h=0$. Thus it removes the coarse-space component of the error and is an orthogonal projection in the energy inner product. With pre- and post-smoothing matrices $S_{\mathrm{pre}},S_{\mathrm{post}}$, the exact two-grid error matrix is $S_{\mathrm{post}}{C}_h{S}_{\mathrm{pre}}$. Effective smoothing on the complementary space and accurate coarse representation of smooth error are both needed for contraction.

A concrete harmonic calculation makes the mechanism quantitative. On the one-dimensional constant-coefficient problem, couple the two harmonics $\theta$ and $\theta +\pi$, put $s={\sin }^2(\theta /2)$ and $c=1-s$, and use the above transfers. Linear interpolation has harmonic weights $(c,s{)}^T$, while the fine operator has diagonal symbol $4\mathrm{diag}(s,c)/h^2$. Therefore its Galerkin correction on this pair isOne pre- and one post-sweep of weighted Jacobi with $\omega =2/3$ gives $E=SCS$, $S=\mathrm{diag}(1-4s/3,1-4c/3)$. This rank-one matrix has eigenvalues zero andThe two-grid Poisson factor with weighted Jacobi is therefore mesh-independent in this model. A periodic constant nullspace must first be removed; a Dirichlet problem has no such nullspace. This is evidence for the design, not a claim that the same factor applies to every PDE or every smoother.

The geometric multigrid V-cycle replaces the exact coarse solve by one recursively defined coarse cycle. A complete implementation on level $\ell$ proceeds as follows:

Boundary values of the solution are imposed on the fine level; the error equation has homogeneous boundary conditions when those values were already imposed exactly. Residuals must be computed after pre-smoothing, not before. The zero initialization on each recursive error solve makes the cycle an error correction rather than an accidental addition of an old coarse solution. For an energy-symmetric cycle usable in conjugate gradient method, choose compatible adjoint transfers and symmetric smoothing, such as forward Gauss-Seidel before correction and its backward sweep afterward.

With bounded stencil size and fixed sweep counts, a sweep, residual calculation or transfer costs $O(N_{\ell })$. Uniform coarsening gives $N_{\ell -1}\approx 2^{-d}{N}_{\ell }$, soIf the full V-cycle has an error reduction factor $q<1$ independent of mesh size, reducing an arbitrary initial algebraic error by a factor $\epsilon$ takes $O(\log (1/\epsilon ))$ cycles and $O(N\log (1/\epsilon ))$ work. Linear work per cycle alone does not prove the mesh-independent factor; it rests on the smoothing and approximation properties just illustrated. Full multigrid starts on a coarse mesh, interpolates its solution to each successive finer mesh, and uses a few cycles per level. When the interpolation accuracy and cycle contraction match the discretization error, this reaches discretization-level accuracy in $O(N)$ work.

Performance depends on the operator. Strong anisotropy can leave error oscillatory in one direction but weakly penalized by the differential operator, defeating point smoothing; line or plane relaxation and selective coarsening address this. Large coefficient jumps demand compatible coarse spaces. Nonrectangular meshes motivate algebraic rather than purely geometric coarsening. Nonlinear problems use a full-approximation formulation instead of the linear residual correction verbatim. Singular Neumann or periodic operators need right-hand-side compatibility and consistent nullspace treatment on every level. Multigrid is fast because its components address different error scales, not because coarsening alone guarantees a good solver.