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

Put $h=1/(m+1)$ and impose homogeneous Dirichlet boundary conditions, as required by the stated continuum eigenvalues. Let $T=\mathrm{tridiag}(1,-2,1)$ be the unscaled one-dimensional Dirichlet discrete Laplacian. For $1\le k\le m$, the vectors $v_j^{(k)}=\sin (jk\pi h)$ satisfyThe identity follows by adding the two neighboring sines; the values at $j=0,m+1$ vanish. These $m$ mutually orthogonal vectors form the discrete sine transform basis. Consequently their tensor products give all $m^2$ two-dimensional modes, not merely a few candidate eigenvectors.

In a compatible grid ordering the axial operator is $A_5=h^{-2}(T\otimes I+I\otimes T)$. The nine-point operator isIndeed, the added Kronecker product supplies the diagonal neighbors with weight $1/6$, subtracts $1/3$ from each axial neighbor and adds $2/3$ to the central coefficient. Thus it has exactly the coefficients of the second discretization. On $v^{(k)}\otimes v^{(l)}$ the eigenvalues areandEquivalently, $h^2{\lambda }_{9,kl}=-10/3+4[\cos (k\pi h)+\cos (l\pi h)]/3+2\cos (k\pi h)\cos (l\pi h)/3$.

For fixed mode numbers $k,l$, expansion of the sines giveswhereas the mixed Kronecker product correction givesThus both eigenvalue approximations are second order for fixed modes. The leading nine-point error includes an additional positive term ${\pi }^4{h}^2{k}^2{l}^2/6$.

There is also an exact comparison, not just an asymptotic one. Since $\sin x<x$ for $0<x<\pi /2$, and both sine factors in the correction are positive,Both estimates lie on the less-negative side of the continuum Laplacian eigenvalue, so the five-point formula is the better eigenvalue estimate for every corresponding interior mode. The five-point versus nine-point Laplacian eigenvalue accuracy comparison does not contradict the useful isotropy or special Poisson accuracy properties of a nine-point stencil. Those properties concern a different error criterion. The fixed-mode expansion is not uniform for modes with $k$ or $l$ comparable to $h^{-1}$; wavelengths close to the mesh scale need not have small relative spectral error.

Work in the standard class of real constant-coefficient linear multistep methods using only first-derivative evaluations. Writewith ${\alpha }_s\ne 0$. The Dahlquist equivalence theorem states that such a method converges on each fixed finite time interval, for sufficiently regular Lipschitz initial value problems and every set of consistent starting values, if and only if it is consistent and zero-stable. Consistency means $\rho (1)=0$ and ${\rho }^{\prime }(1)=\sigma (1)$. Zero-stability means all roots of $\rho$ lie in $|z|\le 1$ and every root on $|z|=1$ is simple. In particular the consistency root at $1$ is simple.

Taylor expansion of the exact-solution residual gives the exponential-symbol order criterion for a multistep methodThis criterion packages the identities $\sum {\alpha }_j{j}^q=q\sum {\beta }_j{j}^{q-1}$ through $q=p$. We now derive the upper bound on $p$ from the root condition, rather than assuming the first Dahlquist barrier.

Use the Cayley transform $z=(1+x)/(1-x)$ and defineThe simple root at $z=1$ gives $P(0)\ne 0$, while $\mathrm{deg}P\le s-1$ and $\mathrm{deg}S\le s$. Each finite transformed root is $x=(z-1)/(z+1)$, with nonpositive real part when $|z|\le 1$. A possible root at $z=-1$ simply lowers the degree of $R$. Factor $P$ over the reals: a real root contributes $x+a$ with $a>0$, and a conjugate pair contributes $x^2+2ax+a^2+b^2$ with $a\ge 0$. Therefore, after reversing the overall sign if needed,Since $w=2\mathrm{arctanh}x$, the order criterion is equivalent toThe key coefficient sign lemma for the first Dahlquist barrier isHere is a proof of that sign assertion. Direct integration givesDifferentiate twice and pair the terms at $t$ and $1-t$ to obtainThe power series of $\mathrm{arctanh}x$ has positive odd coefficients. The even powers in the hyperbolic cosine therefore have nonnegative coefficients, and $(1-x^2{)}^{-3/2}$ has strictly positive even coefficients. Every even coefficient of $F^{\prime \prime }$ is consequently negative. Also $F$ is even and $F(0)=1$, proving the assertion.

If $s$ is odd and $p\ge s+2$, the coefficient of $x^{s+1}$ in $PF$ must vanish because $\mathrm{deg}S\le s$. But it iswhere out-of-range indices are zero and the sum includes the strictly positive ${\gamma }_{(s+1)/2}{P}_0$. This is impossible, so $p\le s+1$. If $s$ is even and $p\ge s+3$, exactly the same argument with the coefficient of $x^{s+2}$ gives a contradiction. Thus $p\le s+2$ in that case. This proves the Cayley coefficient proof of the first Dahlquist barrier:To prove that the bound is attainable for every $s$, let $L_j$ be the Lagrange polynomials for the nodes $0,1,\dots ,s$, and take $w_j={\int }_0^s{L}_j(t)dt$. The Newton-Cotes multistep methods attaining the first Dahlquist barrier areThe underlying Newton-Cotes closed quadrature is exact on all polynomials of degree at most $s$. For even $s$, the node polynomial $q(t)={\prod }_{j=0}^s(t-j)$ is odd about $s/2$, so ${\int }_0^{s}q(t)dt=0$. Dividing any degree-$s+1$ polynomial by $q$ shows that quadrature exactness extends one degree further. Integrating the exact derivative therefore gives method order at least $s+1$ for odd $s$ and at least $s+2$ for even $s$.

Finally $\rho (z)=z^s-1$ has only simple unit-circle roots, so these methods satisfy the root condition. The weights sum to $s={\rho }^{\prime }(1)$, proving consistency. The Dahlquist equivalence theorem gives convergence with consistent starts; the proved upper bound makes the attained formal orders exact; starting values accurate to those orders also give the corresponding global accuracy. Hence the highest convergent order is $2\lfloor (s+2)/2\rfloor$. This result concerns the stated linear first-derivative class, not unrestricted multiderivative or nonlinear formulas.

Let $e$ be the vector of ones. For a general nonautonomous Runge-Kutta method the stages satisfy $k_i=f(t+c_{i}h,y+h{\sum }_j{a}_{ij}{k}_j)$ and the update is $Y=y+h{\sum }_i{b}_i{k}_i$. Expansion at $(t,y)$ givesand thereforeThe exact solution is $y(t+h)=y+hf+h^2(f_t+f_{y}f)/2+O(h^3)$. Matching the independent terms gives the necessary and sufficient second-order conditions with independent Runge-Kutta abscissaeUnder the usual internally consistent Butcher tableau convention $c=Ae$, this reduces to $c=Ae$, $b^{T}e=1$, $b^{T}c=1/2$. If the abscissae are allowed to be independent of the row sums, $c=Ae$ is a sufficient simplifying convention rather than a logically necessary order-two condition.

The original PDF has a sign defect in the lower-left stage coefficient. Put $d=c_2-c_1$. With the printed negative sign, the row sums and the relevant Butcher order conditions areThus the claimed universal second order is false. For example $c_1=0,c_2=1$ gives, on $y^{\prime }=y$, two stages equal to $y$ and the update $Y=(1+h)y$. This is the first-order Forward Euler method. For the literal tableau, order at least two holds exactly when $c_1=1/2$ or $c_2=0$, subject to distinct allowed nodes.

For completeness the literal tableau's A-stability can also be classified algebraically. DefineThe stability function obtained from $R(z)=1+z{b}^T(I-zA{)}^{-1}e$ isIf $D\ge 0$, the denominator has no zero in the closed left half-plane: its reciprocal roots are the stage-matrix eigenvalues, whose positive sum and nonnegative product give positive real parts for nonzero eigenvalues. On the imaginary axis,Consequently for the literal tableau with uniquely solvable stages throughout the left half-plane, the complete parameter criterion isThese explicit inequalities in $c_1,c_2$ are necessary by the small and large imaginary parameters, and sufficient by the maximum modulus principle, including the bound at infinity. For $D=0$ the last inequality forces $r=0$, leaving the stable linear-fractional case. For $D<0$ there is a negative real stage pole. If one defines stability solely through an analytically continued scalar stability function, the only additional possibilities areIndeed these are precisely the cases where the left pole cancels and $R(z)=[1+(1-q)z]/[1-qz]$ is an A-stable theta method. The original stage equations are nevertheless singular at $z=-1/(q-a)$, so they do not give uniquely determined stages there. This distinction prevents a canceled rational factor from concealing an ill-defined stage solve.

The natural intended repair is to replace just the lower-left coefficient by $+c_2^2/(2d)$. This makes the method the collocation Runge-Kutta method on the two nodes. To derive the repaired entries, takeIntegration gives the other three printed entries and the weights unchanged, with the positive lower-left entry. Since $L_1+L_2=1$ and $c_1{L}_1+c_2{L}_2=t$, it follows that $Ae=c$, $b^{T}e=1$ and $b^{T}c=1/2$. Hence the corrected family is always of order at least two.

For the corrected family let $a=(c_1+c_2)/2$, $D=c_1{c}_2/2$ and $r=(1-c_1)(1-c_2)/2$. Its stability function has the same displayed rational form, now with $q=1/2$. The denominator again has no left-half-plane zero, and on the imaginary axis the squared-modulus difference is $(D^2-r^2)y^4$. Because $D,r\ge 0$, the necessary and sufficient inequality is $D\ge r$, orFor sufficiency, $|R|\le 1$ on the imaginary axis and at infinity; apply the maximum modulus principle on growing left half-discs to get the same bound inside. If $D=0$, the condition leaves the endpoint trapezoidal rule and its linear denominator. If $c_1+c_2<1$, every nonzero imaginary argument gives $|R|>1$. This is the two-node collocation A-stability criterion, explicitly distinguished from the literal defective tableau.

Expand the three matrix exponentials and retain terms through degree two. Their product iswhich agrees with $e^{t(A+B)}$ to that order. Keeping the cubic terms gives the more informative defectwhere $[A,B]=AB-BA$ is the commutator. For example, the coefficients of $A^{2}B,ABA,B{A}^2$ in the product are $1/8,1/4,1/8$, instead of the exact $1/6,1/6,1/6$. The cubic defect is generically nonzero. Therefore Strang splitting has order $p=2$. The symmetry $E(-t)=E(t{)}^{-1}$ is consistent with this cancellation of the quadratic defect; for commuting $A,B$ the composition is exact rather than merely second order.

Choose the centered Dirichlet discrete Laplacian in each direction. With $T=\mathrm{tridiag}(1,-2,1)$ put $A=T\otimes I$, $B=I\otimes T$ and $\mu =\Delta t/h^2$. Thus $h^{-2}A$ and $h^{-2}B$ are the physical derivative approximations; the factor $h^{-2}$ has been absorbed into the Courant number in the written update. Equivalently use the scaled derivative matrices and the time parameter $\Delta t$.

For a one-dimensional vector with endpoint values $v_0=v_{m+1}=0$,Thus $A$ and $B$ are real symmetric negative-definite matrices. An orthonormal eigenvector basis shows that all eigenvalues of each exponential lie in $(0,1]$ for $\mu \ge 0$, soIterating gives $\Vert U^n{\Vert }_2\le \Vert U^0{\Vert }_2$ for every nonnegative step parameter, a mesh-independent unconditional stability bound. Multiplying the squared norm by the grid-area weight $h^2$ leaves the same estimate. In this natural tensor-grid choice $AB=BA=T\otimes T$, so $E(\mu ;A,B)=e^{\mu (A+B)}$: the split update is even the exact semidiscrete diffusion flow. Contractivity of a product of these negative symmetric exponentials would prove stability even without commutation.

Write $h=\Delta x$, $k=\Delta t=\mu h$. Every smooth exact solution of the advection equation here is $u(x,t)=F(x+t)$. At $s=x+t$, substitution into the scheme gives the unscaled residualThe constant, linear and quadratic Taylor coefficients vanish. Its cubic coefficient isFor a general smooth function the first-order residual is $2k(u_t-u_x)$, so the proper normalized local truncation error on an exact solution isHence the scheme has formal order two for fixed general positive $\mu$. At $\mu =1/2$ and $\mu =1$, the displayed residual vanishes identically for every $F$, not just to one higher Taylor order: the recurrence transports exactly sampled characteristics when both starting levels are compatible with the exact solution. This formal exactness does not imply stability. At $\mu =0$ there is no positive time step. The two-level recurrence also requires a starting procedure for $U^1$; one initial solution profile alone is not two independent starting levels.

Insert a Fourier mode $U_m^n=G^n{e}^{im\theta }$. With $r=1-2\mu$ and $\eta =e^{i\theta }$, the amplification polynomial isIts two roots have product $-\eta$, of modulus one. Thus if both roots are bounded by one, both must have modulus exactly one. Set $G=e^{i\theta /2}\lambda$. The equation becomesWhen $|r|<1$, both roots have modulus one and their separation obeys $|G_{+}-G_{-}|\ge 2\sqrt{1-r^2}>0$ for every frequency. The companion matrix has bounded entries and can be diagonalized with the two eigenvectors $(G_{+},1{)}^T,(G_{-},1{)}^T$. The uniform separation bounds the inverse eigenvector matrix. Hence every matrix power is bounded uniformly in $n$ and $\theta$, and the Parseval identity proves L2 norm stability for the two starting levels on the whole grid.

If $|r|>1$, frequencies near $\theta =\pi$ give one root outside the unit circle, so the scheme is unstable. If $|r|=1$, the Nyquist polynomial is $(G-r{)}^2$, a repeated unit root. The corresponding companion matrix has a nontrivial Jordan block and solutions containing $n{r}^n$. Thus the endpoints are unstable as well. On the whole line the individual Nyquist plane wave is not square-summable, but Fourier packets localized arbitrarily close to that frequency realize the unbounded power norm by continuity. Therefore the uniform stability of a shifted two-level advection scheme criterion isIn particular $\mu =1/2$ is stably exact with compatible starts, whereas $\mu =1$ is formally exact but unstable under arbitrary perturbations. Checking only the moduli of the roots and retaining both endpoints would miss this essential distinction.

Sparse Gaussian elimination seeks a direct solution while preserving the advantage of a small number of nonzeros. The main issue is that eliminating unknowns creates new interactions between the remaining unknowns. For a symmetric positive-definite system, simultaneous row and column permutation merely changes the elimination order, and Cholesky decomposition or LDL decomposition can proceed without pivoting. Partition one elimination step asThe remaining system has the Schur complement $C-b{b}^T/a$. It is positive-definite because, for any nonzero $z$, minimizing the positive quadratic form of $A$ over the first coordinate gives $z^T(C-b{b}^T/a)z>0$. This both justifies the next positive pivot and identifies exactly where fill can appear.

The matrix graph of a symmetric sparsity pattern has one vertex per unknown and an edge $i-j$ when $a_{ij}$ is structurally nonzero. At pivot $k$, the updatecan change an entry only when both $i$ and $j$ are current neighbors of $k$. Thus the elimination graph operation removes $k$ and makes its remaining neighbors a clique. Edges absent before this update are fill-in. The graph describes structural or generic nonzeros: special cancellations in the numerical values can reduce the actual fill, but cannot justify omitting possible entries in advance.

A useful global description is the fill-path criterion for symmetric elimination. After a set $E$ has been eliminated, remaining vertices $i,j$ are adjacent generically exactly when their original graph contains a path from $i$ to $j$ with all internal vertices in $E$. The assertion starts with the original edges when $E$ is empty. Eliminating the next vertex concatenates two earlier such paths through it; conversely, split any allowed path at the newly eliminated vertex to recover the two earlier connections. This induction explains why sparse long-range couplings may eventually become dense.

The order is decisive. For a star with four leaves, eliminating a leaf only changes the diagonal of its sole neighbor and creates no fill. Repeated leaf elimination preserves this property. Eliminating the central vertex first instead connects all four leaves, creating six edges and a dense four-variable trailing system. This occurs, for example, in a positive-definite star matrix with central diagonal $5$, leaf diagonals $2$ and edge entries $-1$; no cancellation removes the new leaf couplings.

For an arbitrary graph, an ordering creates no fill precisely when each pivot's later neighbors are already a clique, namely a perfect elimination ordering. Such an ordering exists exactly for a chordal graph. A chordless cycle of length at least four illustrates the obstruction: its first removed vertex connects its two previously nonadjacent neighbors. Eliminating in a chosen order adds chords until the completed graph is chordal. Thus sparsity of the original matrix alone does not guarantee sparse factors; the graph and its ordering matter together.

If pivot $j$ has $d_j$ later neighbors, its column contains $d_j+1$ factor entries and its symmetric outer-product update costs $O(d_j^2)$ arithmetic. The total storage is governed by ${\sum }_j(d_j+1)$ and the work by ${\sum }_j(d_j+1{)}^2$. A band ordering with matrix bandwidth $b$ gives the familiar $O(Nb)$ storage and $O(N{b}^2)$ work, but a narrow profile or graph separator can do better than a single global band. The minimum degree algorithm greedily chooses a current low-degree pivot; it often keeps immediate updates small, but it must account for previously created fill and need not minimize total fill.

Nested dissection instead finds a small vertex separator, recursively eliminates disconnected subdomains, and leaves separator variables until last. For a regular two-dimensional grid with $N$ unknowns, a separator has $O(\sqrt{N})$ vertices. Its dense final work is $O(N^{3/2})$, and recursive balanced separation givesAt each level the summed separator-square storage is $O(N)$, giving $O(N\log N)$ factor storage over the logarithmic number of levels. These estimates assume the usual geometric separator structure; they are not guarantees for every sparse graph. Independent subdomains also expose parallel work.

In a practical solver a symbolic factorization phase uses the pattern and ordering to allocate factor storage and schedule updates. For the lower-triangular factor, the elimination tree takes the parent of column $j$ to be its first later row index containing a structural nonzero. Descendant updates flow towards those later pivots, making the tree useful for dependency scheduling and the organization of the factorization. Compressed sparse storage records only entries and indices; the numerical phase then computes the values, followed by forward and backward substitution. The factors can be reused for many right-hand sides.

Finally, structural efficiency must be balanced with numerical stability. For a nonsymmetric system the pattern can be represented by a directed or bipartite graph, and elimination creates edges from incoming to outgoing neighbors. Separate row and column permutations are possible. A fill-reducing column ordering does not ensure a safe pivot: scaling and partial pivoting may be needed to prevent division by a tiny number and large growth. Such row interchanges can change the predicted fill. Positive-definite symmetric problems avoid this conflict; general sparse solvers must manage it explicitly. Graph-based ordering controls factor sparsity, while pivoting controls numerical safety; an effective sparse factorization needs both.

Stability of a numerical method controls how perturbations in initial data, roundoff and local residuals propagate. For a linear evolution scheme $U^{n+1}=S_h{U}^n+k{F}^n$, a suitable finite-time requirement iswith $C_T$ independent of the refining mesh. The norm should approximate the norm of the continuous problem, such as an area-weighted L2 norm. The bound need not imply decay: $C_T=e^{CT}$ is appropriate when the equation itself permits finite-time growth. For a multilevel scheme the state includes all required time levels.

The error recurrence shows why stability matters. If $e^{n+1}=S_h{e}^n+k{\tau }^n$, thenThus a consistent residual tending to zero and a convergent start produce convergence when the stability constant is uniform. Tiny defects need not remain tiny without that bound. The Lax equivalence theorem makes this precise: for a well-posed linear initial value problem and a consistent linear approximation, stability is equivalent to convergence in the compatible solution norms. It is not a theorem for arbitrary nonlinear schemes, and a high formal order does not substitute for stability.

For constant coefficients on a periodic or infinite grid, von Neumann stability analysis diagonalizes spatial translations using Fourier modes. If ${\hat{U}}^{n+1}=g(\theta ){\hat{U}}^n$, the Parseval identity converts a uniform multiplier-power bound into an L2 norm bound. The contraction condition $|g(\theta )|\le 1$ is sufficient; a more general $|g(\theta )|\le 1+Ck$ also yields a finite-time bound. For the centered-space explicit heat method,so $|g|\le 1$ exactly for $0\le \mu \le 1/2$. Larger fixed Courant numbers amplify high frequencies. The Backward Euler method gives instead $g=[1+4\mu {\sin }^2(\theta /2){]}^{-1}$ and is contractive for every $\mu \ge 0$.

For $u_t=u_x$, forward time and centered space give $g=1+i\nu \sin \theta$, with $\nu =k/h$. Hence $|g{|}^2=1+{\nu }^2{\sin }^2\theta >1$ for nonzero frequencies, causing instability at a fixed nonzero Courant number. The directionally correct upwind finite difference scheme uses the neighbor at $m+1$ and hasIt is contractive exactly for $0\le \nu \le 1$. This also has a direct maximum-norm proof: the update is a convex combination of two old values. The continuum domain of dependence must fit inside the numerical one for a local explicit hyperbolic method to converge, giving the Courant–Friedrichs–Lewy condition. That necessary geometric restriction is not by itself a sufficient stability proof.

For systems and multilevel methods, each Fourier multiplier becomes a matrix. All amplification roots must obey the appropriate root bound, but their moduli alone are insufficient. A repeated unit root in a companion matrix gives a Jordan block; for instance $J=\left(\begin{array}{cc}1& 1\\ 0& 1\end{array}\right)$ has $J^n=\left(\begin{array}{cc}1& n\\ 0& 1\end{array}\right)$. Even distinct roots require control of the diagonalizing matrices uniformly in frequency and mesh. The uniform power bound from separated amplification roots supplies such control for a bounded two-root family with a positive gap. Question5 illustrates both sides: strict $0<\mu <1$ gives a uniform gap, whereas the endpoints have a double unit root and are unstable despite unit moduli.

Energy estimates are especially useful when coefficients vary or boundaries prevent Fourier diagonalization. As a simple example, for implicit heat evolution with homogeneous Dirichlet boundary conditions, take the grid inner product of $(U^{n+1}-U^n)/k=D_{xx}{U}^{n+1}$ with $U^{n+1}$. Discrete integration by parts givesThe boundary term vanishes because the boundary values are zero, proving contraction. For variable positive diffusion coefficients the corresponding flux-weighted squared differences remain nonnegative, so the argument extends without a constant-coefficient Fourier symbol. With forcing or lower-order terms, a discrete energy inequality and Gronwall inequality yield a mesh-independent finite-time bound.

A related technique is eigenvalue stability analysis of a finite difference method. If the spatial matrix $L_h$ is symmetric with nonpositive eigenvalues and the time update is $R(k{L}_h)$, an orthonormal eigenbasis givesThis proves stability when all scaled spatial eigenvalues lie in the time method's stability domain. For a nonnormal spatial or amplification matrix, the eigenvalues alone do not control powers: ill-conditioned eigenvectors or Jordan terms can cause transient growth. Direct operator-norm estimates, symmetrizers or energy estimates then provide the missing information. In particular an A-stable scalar time formula does not automatically give mesh-uniform bounds for every ill-conditioned family of spatial matrices.

Boundary conditions and starting procedures are part of the discretization. Interior Fourier analysis may fail to detect a growing boundary mode, and a multilevel method can excite a parasitic mode through a poor start. One must therefore analyze the actual finite-domain operator or include its boundary terms in an energy proof. Stability is a uniform error-propagation bound; Fourier symbols, energy identities, matrix norms and spectral analysis are complementary ways of establishing that bound.