Past exam of the mathematics course of the University of Cambridge / 2015 / iii / Paper 68 / Section A

The characteristic polynomials of a linear multistep method areThe consistency of a numerical method conditions hold for every real $a$:Both roots of a polynomial of $\rho$, namely $1$ and $-1$, are simple and have modulus one. The root condition for a multistep method therefore gives zero-stability for every $a$. The Dahlquist equivalence theorem now gives convergence for every fixed $a\in \mathbb{R}$. As usual, this means convergence of a numerical method on each fixed finite interval for a sufficiently regular ordinary differential equation with a Lipschitz continuous vector field and starting values tending to the exact starting values. When $a\ne 0$, the implicit time-stepping method update is locally uniquely solvable for sufficiently small $h$, for example by a contraction mapping if $|ha|L<1$. Zero-stability does not assert that a large fixed step is suitable for a stiff differential equation.

Expand the exact solution about $t_{n+1}=t$. The unscaled local truncation error isOnly odd powers occur in this centered expansion. Unless $a=1/3$, the first nonzero coefficient is the $h^3$ coefficient, giving order of a numerical method two. For $a=1/3$, that term vanishes but the next coefficient is $-1/90$, giving order of a numerical method four. ThusThese are exact orders for general smooth ordinary differential equations, since the respective first surviving derivatives need not vanish. With starting errors $O(h^p)$, the zero-stability established in part (a) makes the global error $O(h^p)$.

Apply the Dahlquist test equation and put $z=h\lambda$. Every recurrence mode must satisfyThe amplification polynomial of a multistep method, rather than just its root near one, determines absolute stability. Put $B=2a-1$ and use the Cayley transform between the half-plane and diskFor $\zeta \ne -1$, the characteristic equation becomesIf $B>0$, the denominator cannot vanish at a root of the characteristic equation: $1+B{w}^2=0$ would force $w=0$, a contradiction. Hence $\mathrm{Re}z<0$ forces $\mathrm{Re}w<0$, so all amplification roots have modulus less than one. On the imaginary axis the roots have modulus one and are simple: the transformed quadratic has discriminant $4(1-B{z}^2)>0$ for imaginary $z$. At $z=0$ they are the two simple roots $\pm 1$. Also $\zeta =-1$ cannot be a root when $B>0$ and $z\ne 0$, and the leading coefficient $1-az$ cannot vanish in the closed left half-plane.

If $a<1/2$, the root near $-1$ isFor small negative real $z$ it lies below $-1$, violating absolute stability. The endpoint $a=1/2$ deserves separate treatment:One root is always $-1$; the other is the trapezoidal rule multiplier $(1+z/2)/(1-z/2)$. They are distinct for every finite $z$ in the left half-plane. Consequently, with absolute stability understood as the bounded root condition for a multistep method,There is a convention at this reducible endpoint: if A-stability is defined to require every unreduced recurrence mode to decay for $\mathrm{Re}z<0$, the answer is $a>1/2$, since the $(-1{)}^n$ mode persists at $a=1/2$. Canceling the common factor $\zeta +1$ gives the A-stable trapezoidal rule, but cancellation removes an actual starting-error mode of the original two-step recurrence. The fourth-order member $a=1/3$ is outside either A-stability range.

Use the Dirichlet Laplacian eigenvalues and eigenfunctions on $(-1,1)$:An expansion $u={\sum }_j{q}_j(t)e_j$ reduces the wave equation to independent ordinary differential equationsWhen $\alpha <{\lambda }_1$, put ${\omega }_j=\sqrt{{\lambda }_j-\alpha }$; each coefficient isFor finite-energy data $\varphi \in H_0^1(-1,1)$ and $\psi \in L^2(-1,1)$, a uniform spatial bound follows from the energy method, rather than from unjustified absolute summation of the Fourier series. The conserved quantity isThe Poincare inequality gives $\Vert u{\Vert }_2^2\le {\lambda }_1^{-1}\Vert u_x{\Vert }_2^2$. Hence $2E\ge c_{\alpha }\Vert u_x{\Vert }_2^2$ with $c_{\alpha }=1-\alpha /{\lambda }_1>0$ for $0\le \alpha <{\lambda }_1$, and with $c_{\alpha }=1$ for $\alpha <0$. Since the Dirichlet boundary condition gives $u(x,t)={\int }_{-1}^x{u}_x(s,t)ds$, the Cauchy-Schwarz inequality yieldsAt $\alpha ={\lambda }_1$, the first coefficient is $q_1(t)={\varphi }_1+t{\psi }_1$, which is unbounded for some admissible data. For $\alpha >{\lambda }_1$, that coefficient has an exponentially growing component for some admissible data. Thus all-time boundedness of a wave equation with a reaction term requiresAt equality, boundedness for a particular data set requires ${\psi }_1=0$; the remaining modes are bounded by their spectral gap. Above the threshold, every unstable mode must have its growing component canceled, and any zero-frequency mode must have zero initial velocity. There is no condition on $\alpha$ alone for arbitrary specially chosen data.

The printed reference to a limit needs qualification. Bounded oscillations generally have no limit as $t\to \infty$. If actual existence of that limit for every admissible initial datum is required, no real $\alpha$ works: for every $\alpha$, choose $j$ with ${\lambda }_j>\alpha$ and nonzero pure oscillatory data in that eigenfunction. The boxed inequality answers the intended long-time boundedness question.

Set $h_x=1/(N+1)$ and $J=2N+1$. The Dirichlet discrete Laplacian has the eigenvectorsand its negative has the eigenvaluesThe discrete sine transform therefore reduces the method of lines equations to $q_j^{\prime \prime }=(\alpha -{\lambda }_j^h)q_j$. Each coefficient is an oscillator, a linear function, or a hyperbolic function according to the sign of ${\lambda }_j^h-\alpha$. A finite collection of oscillators is bounded in every fixed-grid norm. A zero-frequency mode can grow linearly, and a negative-frequency-square mode can grow exponentially. Consequently,is necessary and sufficient for all-time boundedness of a semidiscrete reaction wave equation on this grid. Equality allows bounded displacement only when the initial velocity has zero projection onto the first eigenvector; above it, the growing components must be absent in every unstable mode. A strictly positive spectral gap also gives a time-independent displacement bound uniform over grids whose gaps are bounded below.

The discrete threshold is strictly smaller than ${\pi }^2/4$, because $\sin s<s$ for $s>0$, and tends to ${\pi }^2/4$ as $N\to \infty$. A coarse grid can therefore introduce long-time growth into a continuum problem with $\alpha <{\pi }^2/4$. This is distinct from displacement stability of a symmetric semidiscrete wave equation on fixed finite time intervals. As in part (a), if existence of the printed limit is required for every initial datum, no $\alpha$ qualifies: any one mode admits oscillatory, secular, or growing initial data, depending on its coefficient.

For the quadratic variational principle for a symmetric positive operator below, use a symmetric operator $L$ on the real Hilbert space $H$. Strict positive definiteness meansA uniformly positive definite symmetric operator satisfies the stronger coercive operator conditionIn this variational setting, “positive definite” is often used for a symmetric operator with this uniform bound. We will state explicitly where the coercive operator bound is needed. For a bounded linear operator defined on all of $H$, symmetry means that the operator is self-adjoint. For an unbounded operator, positivity is imposed on its operator domain, and the variational formulation is made on its form domain.

Symmetry is essential in a real Hilbert space: positivity of the quadratic expression alone does not imply symmetry. For example, $L=I+K$ with nonzero real skew-symmetric matrix $K$ satisfies $⟨Lv,v⟩=\Vert v{\Vert }^2$, but its quadratic functional has derivative involving $I$, not $I+K$. A positive definite symmetric operator supplies both the positivity and symmetry needed in part (b).

First let $L$ be a bounded linear operator on $H$ satisfying the symmetric positive-definiteness convention in part (a). For any direction $w\in H$, expansion of the quadratic functional givesThus vanishing of the first variation in every direction is precisely the weak equation $⟨Lu,w⟩=⟨f,w⟩$ for every $w$. In this whole-space bounded-operator setting it is equivalent to $Lu=f$, the Euler-Lagrange equation.

If $u$ solves that equation, set $v=u+w$. The linear terms cancel:with equality only when $v=u$. Hence the weak solution is the unique global minimizer. Conversely, any minimizer has zero first variation, and therefore solves the weak equation. For a symmetric bounded bilinear form $B$ on a form space $V$, exactly the same calculation gives $I(v)=B(v,v)-2\ell (v)$ and $B(u,w)=\ell (w)$ for every $w\in V$; it does not require an unbounded differential operator to map every $v\in V$ into $H$.

Existence for every $f$ requires an extra hypothesis if “positive definite” means only strict positivity. The coercive operator bound makes the form coercive, so the Lax-Milgram theorem supplies existence and uniqueness. Without that bound, the diagonal operator on sequence space $L(x_n)=(x_n/n)$ on ${\ell }^2$ is symmetric and strictly positive, but $f=(1/n)$ belongs to ${\ell }^2$ and its formal inverse $(1,1,\dots )$ does not. Thus strict positivity alone proves uniqueness and the minimizing property of a solution when one exists, not existence for all $f$.

There is a genuine conflict in the printed coefficient assumptions: no uniformly positive coefficient in a zero-boundary Sobolev space exists. A function $a\in H_0^1(\Omega )$ cannot also obey $a\ge a_{-}>0$ almost everywhere. Indeed, the Lipschitz truncation $g(s)=\min (\max (s,0),a_{-})$ has $g(0)=0$, so Lipschitz truncation preserves zero-boundary Sobolev spaces, giving $g(a)\in H_0^1(\Omega )$. But $g(a)$ would be the nonzero constant $a_{-}$. Its zero gradient contradicts the Poincare inequality in the zero-boundary Sobolev space. Thus the literal coefficient class is empty.

For the meaningful uniformly elliptic problem, take $a\in L^{\infty }(\Omega )$ with $0<a_{-}\le a\le a_{+}$, and impose the Dirichlet boundary condition on the unknown and test functions: $V=H_0^1((0,1{)}^2)$. Additional $H^1$ regularity of $a$ is harmless, but a zero trace for $a$ must be removed. The divergence-form elliptic operator is $L=-\mathrm{div}(a\nabla )$. Integration by parts defines the symmetric bounded bilinear formIn particular,Here $2{\pi }^2$ is the first Dirichlet Laplacian eigenvalue on the unit square; the Poincare inequality follows, for example, by applying the one-dimensional inequality in each coordinate and adding. Thus $B$ is coercive in the gradient norm on $V$, and the Dirichlet realization of $L$ is a positive definite symmetric operator. On its operator domain, $⟨Lv,v⟩=B(v,v)$.

The required functional and weak equation areBecause $f\in L^2$, the right side is a bounded linear functional by the Cauchy-Schwarz inequality and Poincare inequality. The Lax-Milgram theorem gives a unique weak solution, and part (b) proves that it uniquely minimizes $I$. These formulas prove the intended conclusion under the repaired coefficient hypothesis; under the literal hypothesis there is no coefficient to which the conclusion can be applied.

The tableau is the three-stage Lobatto IIIA method. To check its order of a Runge-Kutta method directly, write $e=(1,1,1{)}^T$, $c=Ae=(0,1/2,1{)}^T$, and $C=\mathrm{diag}(c)$. The eight fourth-order conditions for a Runge-Kutta method evaluate toPowers of $c$ here mean componentwise powers. For example, $Ac=(0,1/8,1/2{)}^T$, $A{c}^2=(0,1/24,1/3{)}^T$, and $A^{2}c=(0,1/48,1/6{)}^T$, which make the last three checks immediate. These Butcher order conditions establish order at least four for general nonlinear ordinary differential equations.

It is not order five: part (b)'s stability function has expansionThe fifth coefficient already fails on the Dahlquist test equation. The order is exactly four.

For a Runge-Kutta method, the stability function is $R(z)=1+z{b}^T(I-zA{)}^{-1}e$. Solving the stage equations on the Dahlquist test equation givesThe denominator has roots of a polynomial $3\pm i\sqrt{3}$, both in the open right half-plane. If $z=x+iy$, direct expansion givesSince $|z{|}^2+12>0$, the modulus of $R$ is at most one exactly when $x\le 0$. ThusThe method is A-stable. Its stability function has unit modulus on the imaginary axis and tends to one as $z\to -\infty$, so it is not L-stable.

The criterion for algebraic stability of a Runge-Kutta method is $b_i\ge 0$ and positive semidefiniteness ofThe Runge-Kutta contractivity identity implies that this criterion is sufficient for B-stability, namely contractivity for a dissipative vector field, whenever the stage equations are defined. The implication from algebraic stability to B-stability is the relevant nonlinear theorem; A-stability alone is a linear property.

All weights here are positive, butA positive semidefinite matrix cannot have this negative quadratic value. Hence the method is not algebraically stable. Failure of this sufficient criterion alone would not be a proof of failure of B-stability; no such converse is needed for the question.

For energy stability for variable-coefficient reaction diffusion, let $h_x=1/(M+1)$, let $D_h$ be the Dirichlet discrete Laplacian, and set $V_h=\mathrm{diag}(a_m)$. Use the mesh-weighted inner product $(v,w{)}_h=h_x{\sum }_{m=1}^M{v}_m{w}_m$ with the corresponding discrete L2 norm. For zero boundary values, summation by parts givesFor a solution, or a difference of two solutions, the energy method givesThe Gronwall inequality therefore provesThe constant is independent of $M$, so this is stability of a numerical method on every fixed finite time interval. The bound allows physical growth if $a_{+}>0$; stability here does not mean uniform boundedness as $t\to \infty$.

There is also a uniform maximum-norm bound. At a positive spatial maximum the second difference is nonpositive. Applying this observation to $e^{-a_{+}t}v$ and to its negative gives $\Vert v(t){\Vert }_{\infty }\le e^{a_{+}t}\Vert v(0){\Vert }_{\infty }$. With a perturbing source $r(t)$, the Duhamel principle gives the corresponding initial-data bound plus ${\int }_0^t{e}^{a_{+}(t-s)}\Vert r(s)\Vert ds$ in either norm. Thus the estimate controls accumulated residuals as well as initial perturbations.

The original PDF has $u_m^{n+1}$ on the left of the update; the converted TeX omits its $+1$. Use the Forward Euler method with $k=\Delta t\ge 0$ and $\mu =k/h_x^2$. The explicit time stepping for bounded reaction diffusion update matrix isWhen $0\le \mu \le 1/2$, the heat update $H_h$ has nonnegative stencil weights $\mu ,1-2\mu ,\mu$ and row sums at most one. Hence its induced maximum operator norm is at most one. It is also a symmetric matrix with eigenvaluesall in $[-1,1]$, so its discrete L2 norm operator bound is at most one too. Let $A_{\ast }=\max (|a_0|,|a_{+}|)$. In either of these norms,Iteration gives the mesh-uniform stability estimateThe same estimate applies to differences of solutions. Adding a per-step defect $k{r}^n$ gives $\Vert e^n\Vert \le e^{A_{\ast }T}(\Vert e^0\Vert +k{\sum }_{j<n}\Vert r^j\Vert )$.

This proof only requires $\mu \le 1/2$ and bounded reaction coefficients. It does not assume the full update $E_h$ is nonnegative: a negative $a_m$ can make its middle stencil weight negative when $\mu =1/2$. Nor does it claim contractivity or an all-time bound for every coefficient. The required stability is a bound uniform in the discretization on fixed finite intervals.