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

Take real $\kappa$ and use the L2 norm on the spatial interval. Existence, uniqueness and continuous dependence are the three requirements of Hadamard well-posedness. An energy method supplies the decisive estimate. For a smooth solution with homogeneous Dirichlet boundary conditions, integration by parts givesThe drift contributes only a boundary term, which vanishes. The Poincare inequality $\Vert u_x{\Vert }_2^2\ge {\pi }^2\Vert u{\Vert }_2^2$ further givesApply the same argument to the difference of two solutions to obtain uniqueness and continuous dependence on the initial data.

For existence, use the Dirichlet gauge transform for constant drift: $w=e^{\kappa x/2}u$ satisfies $w_t=w_{xx}-{\kappa }^{2}w/4$ with zero boundary values. Expanding $w_0$ in its Fourier sine series givesFor $u_0\in L^2(0,1)$ the series defines a solution continuous in $L^2$ down to $t=0$ and smooth for positive time; multiplication by the fixed bounded exponentials preserves this interpretation. Its energy estimate follows by approximation with smooth initial data. For a classical solution at the initial corners, require the usual smoothness and boundary compatibility instead. The problem is well posed in $L^2$, with a contraction estimate independent of the initial data.

Let $d=\Delta x$, impose $u_0=u_{M+1}=0$, and write the method of lines system as $U^{\prime }=LU$, whereThus $D_2$ is a negative definite symmetric matrix and $D_1$ is a skew-symmetric matrix. Use the mesh-weighted Euclidean norm $\Vert U{\Vert }_d^2=d{\sum }_{m=1}^M|u_m{|}^2$. Discrete summation by parts yields the centered Dirichlet drift-diffusion energy identityConsequentlyThe same estimate controls perturbations and is uniform in the number of grid points and in the fixed drift coefficient. Finite-dimensional linear ODE theory guarantees existence, so this proves stability of a numerical method for the semidiscretization.

The factor $d^{1/2}$ simply rescales the vector norm and does not change the induced matrix norm. Equivalently the symmetric part is $(L+L^{\ast })/2=D_2$, whose largest eigenvalue is $-4d^{-2}{\sin }^2[\pi /(2(M+1))]<0$. This is the Euclidean logarithmic norm, rather than generally the spectral abscissa of a nonnormal matrix. No periodic Fourier mode assumption has been made: the zero endpoint terms are part of the proof. In particular positivity of both off-diagonal coefficients is not needed for this $L^2$ stability result.

Put $k=\Delta t>0$ and retain the same matrix $L$. The Crank--Nicolson method isTake the real mesh-weighted inner product with $U^{n+1}+U^n$. The cross terms cancel, and the identity from part (b) givesThe implicit system is uniquely solvable: if $(I-kL/2)V=0$, thenso $V=0$. Therefore its dissipative Cayley-transform contraction satisfiesIterating proves unconditional stability for every $\mu =k/d^2>0$; the perturbation bound is one and does not depend on the mesh or time-step ratio.

The printed hint's exponential estimate is valid with the logarithmic norm, but its proposed bound by $r(k\alpha [L])$ is not a general inheritance principle. For the trapezoidal stability function $r(z)=(1+z/2)/(1-z/2)$, that real number can even be negative when $k\alpha [L]<-2$, whereas a norm is nonnegative. The direct energy proof above establishes the required result without that assertion.

The characteristic polynomials of a linear multistep method areTo determine formal order, substitute a smooth exact solution and expand about the first time level. The exponential-symbol order criterion for a multistep method collects precisely the same coefficients:The constant, linear and quadratic coefficients vanish for every $a$. The cubic coefficient vanishes only at $a=-1/5$, where the quartic coefficient is $-1/30\ne 0$. HenceHere order means the exact-solution step residual is $O(h^{p+1})$. It is a formal consistency result, not a convergence assertion. In particular at $a=1$ both ${\rho }^{\prime }(1)$ and $\sigma (1)$ vanish, and the double root at one destroys zero-stability; cancelling its common factor gives a different, first-order recurrence with an additional integration constant left unspecified by the original formula.

The Dahlquist equivalence theorem states that a consistent linear multistep method is convergent for suitably consistent starting values exactly when it is zero-stable. The root condition for a multistep method requires every root of $\rho$ to lie in the closed unit disk, with every unit-modulus root simple.

The roots are $1$ and $a$. Thus $|a|\le 1$ is necessary; $a=1$ is excluded because it gives a double root at one. At $a=-1$, the two unit roots are distinct, so the endpoint is allowed. Combined with part (a), this provesAssume a locally Lipschitz vector field, a smooth solution on the fixed time interval, a nearby solvable implicit branch and starting errors of the required order. The global order is three at $a=-1/5$ and two at the other convergent parameter values. Outside this interval, zero-step perturbations already grow through either an exterior root or a unit-root polynomial factor.

Use the standard A-stability convention including the root condition at $z=h\lambda =0$. The amplification polynomial of a multistep method isZero-stability first restricts the possible parameters to $-1\le a<1$. For $-1<a<0$, as $z$ tends to negative infinity, one amplification root tends to $(3a-1)/(1+a)$, whose modulus exceeds one. For $a=-1$, the characteristic equation is ${\zeta }^2-2z\zeta -1=0$; large negative real $z$ likewise gives an exterior root. Thus $a\ge 0$ is necessary.

For sufficiency take $0\le a<1$. The leading coefficient cannot vanish in the closed left half-plane. On the unit circle $\zeta =e^{i\theta }$ the boundary-locus test for multistep A-stability uses $z=\rho (\zeta )/\sigma (\zeta )$ and giveswhere the denominator is nonzero. When $a=0$, the exceptional zero of $\sigma$ at $\zeta =-1$ is not a zero of $\rho$ and therefore cannot be an amplification root for a finite $z$. Near $z=0$ on the negative real axis, the root issuing from one is $1+z+O(z^2)$ and is strictly inside the disk; the other root is close to $a$ and is also inside. Roots vary continuously, cannot escape through infinity because the leading coefficient is nonzero, and cannot cross the unit circle anywhere in the open left half-plane by the displayed boundary formula. Thus every root remains inside there. Continuity gives the boundary case; for $a>0$ a unit root on the imaginary axis can occur only at $z=0$, where it is simple. For $a=0$, cancellation of the harmless zero root leaves the trapezoidal rule, whose amplification factor has modulus at most one.

ThereforeThe third-order member is convergent but not A-stable, consistent with the Second Dahlquist barrier. At $a=1$, $P_z=(\zeta -1)[(1-z)\zeta -1]$ has a permanent unit root and a double root at $z=0$: the unreduced method is not zero-stable and is not A-stable under the stated convention. Testing only the open half-plane while overlooking its zero-step behavior would give a weaker conclusion.

Let $p$ be the collocation polynomial of degree at most $s$ over one step, with $p(t_n)=y_n$, and let $Y_i=p(t_n+c_{i}h)$. Its derivative, of degree at most $s-1$, is fixed by its values at the distinct nodes. Lagrange interpolation therefore givesIntegrating from zero to $\tau$ yieldsAt $\tau =c_i$ this gives the stage equations of an implicit Runge-Kutta method; at $\tau =1$ it gives the update. ConsequentlyConversely, stages satisfying these equations define the integrated polynomial displayed above. It takes the stage values at the nodes and has the required derivative there, so it solves the collocation equations. This proves equivalence for every common solution branch, not merely equality on the scalar test equation. Also ${\sum }_j{a}_{ij}=c_i$, since the Lagrange polynomials sum to one. Stage existence or uniqueness requires the usual implicit-solvability assumptions; for a Lipschitz vector field a sufficiently small step gives a contraction. This is the collocation Runge-Kutta method construction.

The Lagrange interpolation polynomials for these nodes areTheir integrals give the Lobatto IIIA method tableauTo verify its nonlinear order, not just its scalar linear order, set $e=(1,1,1{)}^T$, $c=Ae$ and $C=\mathrm{diag}(c)$. Direct multiplication verifies the fourth-order conditions for a Runge-Kutta method:Powers of $c$ are componentwise. These eight Butcher order conditions establish order at least four for general smooth ODEs. On the Dahlquist test equation, elimination of the stages givesA nonzero fifth-order step defect rules out order five. Thus the method has exactly order four.

Write $P(z)=1+z/2+z^2/12$ and $Q(z)=1-z/2+z^2/12$. The zeros of $Q$ are $3\pm i\sqrt{3}$, so the stability function has no pole in the closed left half-plane. A direct modulus calculation givesFor $\mathrm{Re}z\le 0$ this is nonnegative, and hence $|R(z)|=|P/Q|\le 1$. Therefore the Lobatto IIIA method is A-stable. It is not L-stable, because $R(z)\to 1$ for large negative real $z$; unconditional scalar stability need not strongly damp the stiffest modes.

By orthogonal diagonalization of a real symmetric matrix, write $A=Q\mathrm{diag}({\lambda }_1,\dots ,{\lambda }_d)Q^T$, with orthogonal $Q$. Its matrix exponential has the same eigenvectors and positive eigenvalues $e^{t{\lambda }_j}$. Orthogonal invariance of the induced Euclidean norm gives the exact identityThis proves the requested inequality with equality. If a real number $\beta$ gave the bound for every $t\ge 0$, evaluating on a unit eigenvector for ${\lambda }_{\max }(A)$ at any $t>0$ would give $e^{t{\lambda }_{\max }}\le e^{t\beta }$, so $\beta \ge {\lambda }_{\max }$. Thus the stated exponent is the smallest possible. For a symmetric matrix the spectral abscissa and Euclidean logarithmic norm coincide, unlike the general nonsymmetric case in Question 1.

Set $S=A+B$ and $E(t)=F(t)-e^{tS}$. Since $F(0)=I$, $E(0)=0$. Differentiate the ordered exponential products, using that each matrix commutes with its own exponential:Here the commutator convention is $[X,Y]=XY-YX$; this fixes both signs. The error satisfies $E^{\prime }=SE+D$, so the variation-of-constants formula givesThis is the symmetrized exponential-splitting defect identity. Symmetry of $A,B$ was not needed for the identity itself; it will be used to bound their exponentials in part (c).

Let $\beta =\mu [A]+\mu [B]$ and $\gamma =\mu [A+B]$, where all three matrices are symmetric. Submultiplicativity and the triangle inequality giveand similarly the other commutator term is bounded by $2\Vert B{\Vert }_2{e}^{x\beta }$. Apply part (a) also to $A+B$ in the integral from part (b). The outer factor one-half cancels these twos, leavingThe exponential divided difference evaluates this integral. Thuswhen $\beta \ne \gamma$, andwhen they coincide. The second expression is both the direct equal-exponent integral and the continuous limit of the first. The Rayleigh-Ritz variational principle also gives $\gamma \le \beta$, although the integral computation does not require a strict inequality. These are valid coarse norm bounds; the cancellation between the two products can make the actual small-step error substantially smaller.

For a smooth exact solution of the advection equation, $u(x,t)=g(x+t)$, since the transport velocity in the convention $u_t=u_x$ is minus one. Put $d=\Delta x$, $k=\mu d$ and move the scheme's right side to the left. Substitution gives the exact-solution residualIts constant, linear and quadratic Taylor expansion terms cancel. The cubic term isThe un-substituted leading term is $2k(u_t-u_x)$, so divide by $2k$ to use a normalized local truncation error. For a fixed positive Courant ratio,Thus the generic method is second order, provided stability and a compatible second-order starter are supplied.

There are two special ratios rather than an unnoticed higher generic order. At $\mu =1/2$, the scheme is $u_m^{n+1}=u_{m+1}^{n-1}$; both sides lie on the same exact characteristic, and the residual vanishes identically. At $\mu =1$ the previous-level term cancels the unshifted current term for exact data, again producing zero residual on every exact characteristic solution. The first special case is stable, while the second is not stable for arbitrary two-level perturbations, as part (b) shows. Exact propagation of specially initialized data is not a substitute for stability.

Use the Fourier transform for the spatial Cauchy problem, or its periodic analogue, and regard the two starting levels as independently perturbed data. A Fourier mode with spatial factor $e^{im\theta }$ has amplification roots $G$ satisfyingPut $c=2\mu -1$ and $G=e^{i\theta /2}q$. ThenIf $|c|<1$, both roots have modulus one and their separation is bounded below uniformly in frequency:The uniform power bound from separated amplification roots now controls the two-level companion matrix for every time step. Its entries are uniformly bounded, and its eigenvector conditioning is bounded by the reciprocal root gap. The Parseval identity transfers this frequency-uniform bound to the spatial ${\ell }^2$ norm. This proves stability, rather than merely checking each root's modulus.

If $|c|>1$, the frequency $\theta =\pi$ has a root outside the unit disk, so there is exponential instability. If $|c|=1$, the two roots at $\theta =\pi$ coincide on the unit circle. The companion matrix is not a scalar matrix and has a nontrivial Jordan block; its powers grow linearly in the number of steps. Frequencies arbitrarily near that value produce the same lack of a uniform bound for localized Fourier packets, so this also invalidates Cauchy ${\ell }^2$ stability, not only periodic plane-wave stability. At $\mu =1$ the double amplification root is $-1$, and at $\mu =0$ it is $1$.

Therefore the full two-level stability range for a fixed positive Courant ratio isThe endpoint $\mu =0$ is moreover not a positive time step. Bounds deteriorate as $\mu$ approaches either endpoint; the displayed range is not a uniform claim over ratios arbitrarily close to one. A prescribed starter that removes one special parasitic component can change behavior for selected initial data, but does not establish the requested unrestricted two-level stability.

A stiff differential equation can contain rapidly decaying modes as well as a slowly changing solution of interest. An explicit scheme may need a very small step solely to keep those already small fast modes from numerical growth. A-stability addresses this restriction through the Dahlquist test equation $y^{\prime }=\lambda y$, whose exact solution decays when $\mathrm{Re}\lambda <0$.

For a one-step method, write $y_{n+1}=R(z)y_n$, $z=h\lambda$. Its linear stability domain consists of the points where the update is defined and $|R(z)|\le 1$. A-stability means that this domain contains the closed left half-plane. This is a scalar linear stability property, not an order condition, an existence theorem for an implicit solve, or a general nonlinear contractivity assertion.

For a rational approximation to the exponential, consistency requires $R(z)=1+z+O(z^2)$, and order $p$ requires $R(z)-e^z=O(z^{p+1})$. Poles must be absent from the relevant half-plane. The Forward Euler method has $R=1+z$ and the disk $|1+z|\le 1$, so negative real modes require $h|\lambda |\le 2$; it is not A-stable. The Backward Euler method has $R=(1-z{)}^{-1}$, and $|1-z|\ge 1$ on the left half-plane, so it is A-stable. The trapezoidal rule and the implicit midpoint rule both have, on the scalar linear problem,They are second-order and A-stable despite being different nonlinear methods. No nonconstant polynomial approximation is A-stable, because its modulus grows without bound on the negative real axis. In particular every nontrivial explicit Runge-Kutta method has a bounded-step stability restriction.

A-stability prevents growth but need not eliminate extremely stiff modes: the trapezoidal factor tends to minus one as $z\to -\infty$. L-stability additionally requires $R(z)\to 0$ within the left half-plane. Backward Euler is L-stable; trapezoidal and implicit midpoint are not. The distinction explains persistent alternating transients in otherwise stable Crank-Nicolson calculations. Accuracy still constrains useful step sizes even for an A-stable method.

A linear multistep method is characterized by polynomials $\rho ,\sigma$, and its scalar modes satisfy $\rho (\zeta )-z\sigma (\zeta )=0$. There is generally no single scalar amplification factor: every root must lie in the closed unit disk and unit roots must be simple. At $z=0$, this is the zero-stability root condition. Consistency plus zero-stability gives convergence by the Dahlquist equivalence theorem, with suitable starting values. A-stability demands the amplification root condition throughout the left half-plane, including the zero-step condition. The Second Dahlquist barrier says that an irreducible A-stable linear multistep method has order at most two, and an explicit linear multistep method cannot be A-stable.

Examples include backward Euler, trapezoidal, and the second-order backward differentiation formulaIts zero-step roots are $1,1/3$. Its unit-circle boundary quotient has real part $(1-\cos \theta {)}^2\ge 0$, and the leading coefficient is nonzero on the left half-plane, giving A-stability by the same continuation argument as in Question 2. Conversely the third-order member of that question fails A-stability. A formal high order without zero-stability, or after silently cancelling a problematic unit factor, does not evade the barrier.

For an implicit Runge-Kutta method, the stages on the test equation satisfy $Y=e{y}_n+zAY$. Whenever these stages are solvable, its stability function isThis rational function is used to test scalar A-stability; the Butcher order conditions separately establish nonlinear order. The Lobatto IIIA method in Question 3 is an A-stable fourth-order example with $R$ the diagonal degree-two Padé approximant. More generally the $s$-stage Gauss--Legendre Runge-Kutta method has order $2s$ and is A-stable, with the diagonal Padé approximation to $e^z$; its stiff-limit factor is nonzero. The $s$-stage Radau IIA method has order $2s-1$ and is L-stable. Thus Runge-Kutta stages permit arbitrarily high A-stable order, unlike irreducible linear multistep formulas.

For a normal matrix, scalar amplification bounds transfer through unitary diagonalization. For a non-normal matrix, eigenvalue information alone does not assert Euclidean contraction: transient growth and eigenvector conditioning can matter. Direct energy method arguments, such as the dissipative Cayley-transform proof in Question 1, use the symmetric part and are more informative in that setting.

For nonlinear dissipative vector fields, the stronger B-stability property asks that distances between two numerical solutions not increase. A standard sufficient condition for Runge-Kutta methods is algebraic stability: $b_i\ge 0$ andis positive semidefinite. It follows from the Runge-Kutta contractivity identity when the implicit stages exist. The three-stage Lobatto IIIA method is A-stable but has the first diagonal entry of that matrix equal to $-1/36$, so scalar A-stability should not be confused with that sufficient nonlinear condition. Stage solvability of an implicit Runge-Kutta method remains a separate issue and computational cost of the implicit equations must be considered. A-stability removes a scalar decay-mode stability restriction; order, stiff damping, nonlinear stability and solvability remain distinct properties.

The finite element method replaces an infinite-dimensional variational problem by a finite-dimensional one, usually using functions that are polynomial on small mesh elements. Ritz method and Galerkin method describe how the discrete equations are selected; neither intrinsically requires piecewise polynomials, though those spaces make assembly local and sparse.

For a concrete realization of the two-point problem, choose homogeneous Dirichlet boundary conditions at both endpoints of $(0,1)$. The differential expression alone is not a complete boundary value problem, so this boundary choice is an explicit assumption. Let $p\in W^{1,\infty }$ with $p\ge p_0>0$, $q\in L^{\infty }$ with $q\ge 0$, and $f\in L^2$. Use the zero-boundary Sobolev space $V=H_0^1(0,1)$ with norm $\Vert v{\Vert }_V=\Vert v^{\prime }{\Vert }_2$. By the Poincare inequality this is equivalent to its full first-order Sobolev norm. Multiplying by a test function and applying integration by parts givesThe endpoint terms vanish by the essential boundary condition. The form is a symmetric bilinear form, bounded byand it is a coercive bilinear form:The load is a bounded linear functional. The Lax-Milgram theorem gives a unique weak solution and a bound $\Vert u{\Vert }_V\le \Vert \ell {\Vert }_{V^{\ast }}/p_0$. Under the stated coefficient regularity it solves the differential equation in the usual weak sense.

The Ritz method minimizes the energyover $V$, or over a conforming finite-dimensional subspace $V_h\subset V$. Its first variation is $a(v,w)-\ell (w)$; symmetry and coercivity make the functional strictly convex. In fact, if $u$ is the weak solution,Thus the energy minimizer is exactly the weak solution. The Galerkin method instead directly asks for $u_h\in V_h$ with $a(u_h,v_h)=\ell (v_h)$ for all $v_h\in V_h$. For this symmetric coercive problem, Ritz-Galerkin equivalence for a symmetric coercive form shows that the discrete minimizer and Galerkin solution coincide. For a nonsymmetric problem, a Galerkin formulation can still apply but the same quadratic minimization generally cannot represent its full bilinear form; appropriate coercivity or inf-sup conditions must then justify the chosen spaces.

Subtract the continuous and discrete weak equations to obtain Galerkin orthogonality $a(u-u_h,v_h)=0$. Since $a$ defines the energy norm $\Vert v{\Vert }_a=\sqrt{a(v,v)}$, the Pythagorean identity givesIn a general reference norm, the Céa lemma yieldsThe proof uses $a(e,e)=a(e,u-v_h)$ together with coercivity and continuity. Dense approximation spaces therefore imply convergence; the estimate reduces a numerical error problem to an approximation problem.

For implementation, partition the interval by nodes $0=x_0<\cdots <x_N=1$. Choose continuous piecewise-linear functions and the interior piecewise-linear hat functions ${\varphi }_j$ as a basis, with the two boundary coefficients prescribed. Write $u_h={\sum }_j{U}_j{\varphi }_j$. The coefficients satisfyThe stiffness matrix is symmetric positive definite: $U^{T}KU=a(u_h,u_h)>0$ for nonzero interior coefficients. The local supports make it tridiagonal in this one-dimensional linear-element case.

On an element of length $h_e$, for constant element coefficients $p_e,q_e$, the two-node matrix and constant-load vector areFor variable coefficients, integrate $p,q,f$ against the same shape functions rather than silently replacing them by constants. Element contributions are added at shared nodes, and boundary values are eliminated or lifted into the right side. A sparse direct solve or a suitable iterative method gives the nodal coefficients. Numerical quadrature should preserve the required accuracy and, with positive coefficients and weights, the positive energy structure.

A finite element interpolation estimate gives ${\inf }_{v_h}\Vert u-v_h{\Vert }_{H^1}\le Ch\Vert u{\Vert }_{H^2}$ for these elements, where $h$ is the largest element size and the solution has the stated regularity. The Céa lemma then gives first-order energy error. If the dual elliptic problem has $H^2$ regularity, the Aubin–Nitsche duality argument gives second-order $L^2$ error: for $e=u-u_h$, solve $a(v,z)=(e,v{)}_{L^2}$, then use $a(e,z)=a(e,z-I_{h}z)$ to gain one further factor of $h$. Higher-order elements improve rates when the solution is sufficiently smooth.

Nonzero Dirichlet data are handled by a boundary lifting and an affine trial space. Neumann boundary conditions enter naturally through the integrated boundary term; Robin terms modify both the form and the load. With pure Neumann conditions and $q=0$, constants lie in the kernel: existence needs the load/flux compatibility condition, uniqueness needs a mean constraint, and the preceding Dirichlet coercivity argument cannot simply be reused. Conforming approximation, boundary conditions and coercivity are the ingredients connecting the finite-element construction to a justified Ritz or Galerkin solution.