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

Apply the Runge-Kutta method to the Dahlquist test equation and put $z=h\lambda$. Its stability function isSince $A$ is an invertible matrix, $Q$ has polynomial degree $s$. The adjugate matrix formula gives $\mathrm{deg}P\le s$, but the additional condition givesso actually $\mathrm{deg}P\le s-1$. The specified order of a numerical method implies $R(z)=e^z+O(z^{2s})$ near zero. Therefore the Padé identification of a Radau stability function makes $R$ the $[s-1/s]$ Padé approximant to the exponential function: its numerator and denominator have the required degrees, and its Taylor series agrees through degree $2s-1$.

In the permitted A-stability of near-diagonal exponential Padé approximants, take $m=s-1$, $n=s$. Then $n-2\le m\le n$, so the scheme is A-stable. The same stability function also tends to zero at infinity, giving L-stability.

For one stage, first-order consistency of a numerical method gives $b_1=1$. The condition $b^T{A}^{-1}1=1$ then gives $a_{11}=1$. Hence the scheme is the Backward Euler method:It is a one-node collocation Runge-Kutta method with $c_1=1$. Specifically, choose a linear polynomial $p(t)$ on $[t_n,t_n+h]$ with $p(t_n)=y_n$ and impose $p^{\prime }(t_n+h)=f(t_n+h,p(t_n+h))$. Its derivative is constant, so $p(t_n+h)-p(t_n)=h{p}^{\prime }(t_n+h)$, exactly the Backward Euler method update.

Write $h=\Delta x$, $k=\Delta t$, and assume $a,u$ are sufficiently smooth for the following Taylor expansions. The symmetric half-grid diffusion consistency calculation giveswhereThe normalized local truncation error of the full update is consequentlyThe PDF prints an advection equation, $u_t=a(x)u_x$, whereas this stencil approximates conservative diffusion, $u_t=(a{u}_x{)}_x$. For the literal printed equation the leading defect is $(a-a^{\prime })u_x-a{u}_{xx}$, which generally does not vanish: the method is inconsistent, with an $O(1)$ normalized defect.

For the intended variable-coefficient conservative diffusion equation, the first two terms cancel and ${\mathcal{T}}_{h,k}=O(k+h^2)$. Under the parabolic step scaling $k=O(h^2)$, this is $O(h^2)$, while the unnormalized one-step local truncation error is $k{\mathcal{T}}_{h,k}=O(h^4)$. These two orders use different normalizations. With numerical stability and compatible smooth initial and boundary conditions, the intended diffusion problem has global error $O(k+h^2)$ on a fixed time interval. Positivity and boundedness of $a$ alone do not justify the displayed Taylor expansion; coefficient regularity is also needed for that error formula.

Put $r=k/h^2$. The update isUnder $k\le h^2/(2\beta )$, all three weights are nonnegative and sum to one. With the imposed zero Dirichlet boundary conditions, this is a convex combination including any boundary values. Thus the monotone half-grid diffusion update givesDifferences of two computed solutions obey the same estimate, proving numerical stability in the discrete maximum norm, uniformly in $h,k,n$.

There is also an L2 norm proof. The diffusion matrix $L_h$ is symmetric and satisfiesSince $(v_{m+1}-v_m{)}^2\le 2(v_{m+1}^2+v_m^2)$, every eigenvalue lies in $[-4\beta /h^2,0]$. Hence the eigenvalues of $I+k{L}_h$ lie in $[-1,1]$ under the same step restriction, and the spectral theorem for real symmetric matrices gives an L2 norm contraction. These stability estimates are valid for the printed recurrence, although they cannot repair its inconsistency with the printed advection equation.

Expanding the symmetric spatial stencil by the Taylor theorem givesConsistency of a numerical method and cancellation of the $h^2$ term require $\alpha +2\beta +2\gamma =0$, $\beta +4\gamma =1$, and $\beta +16\gamma =0$. Solving this linear system givesThe resulting fourth-order centered second derivative hasThe nonzero next coefficient proves that the highest spatial order of a numerical method is four for this three-parameter stencil.

For the method of lines in part (a), the Fourier symbol of the spatial operator isEach Fourier mode therefore evolves by $\hat{U}(\theta ,t)=e^{t{\lambda }_h(\theta )}\hat{U}(\theta ,0)$, whose multiplier has modulus at most one. The discrete Parseval identity givesThus the fourth-order centered second derivative gives an unconditionally stable semidiscrete evolution in the L2 norm. No time integration method has yet been chosen, so no time-step restriction belongs to this conclusion.

An arbitrary L2 function does not have well-defined point samples. For the stated initial data one may use cell-average projection, $U_m(0)=h^{-1}{\int }_{x_m-h/2}^{x_m+h/2}{u}_0(x)dx$. The Cauchy-Schwarz inequality gives $h{\sum }_m|U_m(0){|}^2\le \Vert u_0{\Vert }_{L^2}^2$, providing a bounded initialization for the stability estimate.

Let $V$ be a real Hilbert space, $a:V\times V\to \mathbb{R}$ a bounded bilinear form, and $\ell \in V^{\prime }$ a bounded linear functional. Assume constants $C<\infty$, $\alpha >0$ satisfyThe second hypothesis says that $a$ is a coercive bilinear form. The Lax-Milgram theorem then supplies a unique weak solution $u\in V$ with $a(u,v)=\ell (v)$ for every $v\in V$, and $\Vert u{\Vert }_V\le \Vert \ell {\Vert }_{V^{\prime }}/\alpha$.

For a finite-dimensional trial linear subspace $V_h\subset V$, the same boundedness and coercivity constants apply to the restricted form. Since a finite-dimensional subspace is closed, $V_h$ is itself a Hilbert space. Thus there is a unique Galerkin method solution $u_h\in V_h$ satisfying $a(u_h,v_h)=\ell (v_h)$ for every $v_h\in V_h$. The Galerkin orthogonality $a(u-u_h,v_h)=0$ gives the Céa lemma:In particular, approximation by the trial spaces implies numerical convergence. For a complex Hilbert space, replace bilinearity by a sesquilinear form, take $\ell$ antilinear in the test variable, and require $\mathrm{Re}a(v,v)\ge \alpha \Vert v{\Vert }_V^2$.

Use the zero-boundary Sobolev space $V=H_0^1(0,1)$ with $\Vert v{\Vert }_V=\Vert v^{\prime }{\Vert }_{L^2}$. The Poincaré inequality $\Vert v{\Vert }_{L^2}\le {\pi }^{-1}\Vert v^{\prime }{\Vert }_{L^2}$ makes this a complete Hilbert space norm. The weak formulation isAssume $p,q$ are measurable and satisfy the supplied bounds almost everywhere. The Cauchy-Schwarz inequality and Poincaré inequality giveBecause $p\ge p_0>0$ and $q\ge 0$,This is the coercive variable-coefficient Dirichlet form. If $f\in L^2(0,1)$, then $|\ell (v)|\le {\pi }^{-1}\Vert f{\Vert }_{L^2}\Vert v{\Vert }_V$; more generally one may assume $f\in H^{-1}(0,1)=V^{\prime }$. Some such source regularity is necessary to invoke the theorem. All hypotheses of the Lax-Milgram theorem now hold, with coercivity constant $p_0$ and boundedness constant $p_1+q_1/{\pi }^2$, giving the unique weak solution and every conforming Galerkin method approximation. This argument does not require a pointwise derivative of $p$.

Along the exact ordinary differential equation solution, the chain rule gives $g(y)=f^{\prime }(y)f(y)=y^{\prime \prime }$. Insert the exact solution into the symmetric two-step two-derivative formula and expand about the middle time $t=t_{n+1}$. Its unnormalized local truncation error isThe coefficients of $h{y}^{\prime }$, $h^3{y}^{\prime \prime \prime }$ and $h^5{y}^{(5)}$ are respectivelyThe centered residual is an odd function of $h$, so every even power cancels. The next coefficient isThus the method has order six. Its zero-stability roots at $h=0$ are the simple roots $1,-1$. With sufficiently accurate starts and smooth derivative evaluation maps, convergence of a zero-stable multiderivative method gives sixth-order error on fixed time intervals. The first Dahlquist barrier for ordinary linear multistep methods does not apply, because this is a multiderivative multistep method.

The stability convention matters. The Cambridge numerical analysis notes define a strict linear stability domain through decay $y_n\to 0$. For a multistep recurrence, robust decay requires every amplification root to lie strictly inside the unit disk. With that convention we can prove the stronger result that this method's strict linear stability domain is empty, and therefore bounded. Mere boundedness of oscillatory sequences gives a different answer, described below.

On the Dahlquist test equation, $g(y)={\lambda }^{2}y$. With $z=h\lambda$, the amplification polynomial of a multistep method isExclude $D=0$, where the intended two-step update cannot be solved. If both polynomial roots have modulus below one, their product requires $|A|<|D|$. Butso necessarily $\mathrm{Re}z<0$. The complex quadratic Schur criterion also requiresIts left-hand side equals $32|\mathrm{Re}z|(|z{|}^2+15)$, while for $\mathrm{Re}z<0$ its right-hand side is only $28|\mathrm{Re}z|(|z{|}^2+15)$. The inequality is impossible. ThusFor comparison, the bounded-root stability set of the symmetric two-derivative formula is unbounded. On $z=iy$, rotate the amplification polynomial into a real reciprocal quadratic. Its two roots are distinct and on the unit circle precisely whenequivalently $y^2<(45-15\sqrt{5})/2$ or $y^2>(45+15\sqrt{5})/2$. In particular, arbitrarily large imaginary parameters give bounded, undamped numerical oscillations. At the two equality thresholds, repeated unit-circle roots produce unbounded solutions. Hence the question's bounded-domain claim holds for the strict decay convention; it does not hold if the phrase is interpreted using the non-growing root condition for a multistep method.

A constant-step linear multistep method has the formwith ${\alpha }_k\ne 0$. The characteristic polynomials of a linear multistep method encode both accuracy and error propagation. The method has order of a numerical method $p$ when insertion of a sufficiently smooth exact solution leaves an unnormalized local truncation error $O(h^{p+1})$. The exponential-symbol order criterion for a multistep method expresses this as $\rho (e^z)-z\sigma (e^z)=O(z^{p+1})$. In particular, first-order consistency of a numerical method requires $\rho (1)=0$ and ${\rho }^{\prime }(1)=\sigma (1)$.

Accuracy alone is insufficient. For $f=0$, perturbations obey the recurrence $\rho (E)e_n=0$. The root condition for a multistep method requires every root of a polynomial $\rho$ to lie in the closed unit disk, with each unit-circle root simple. This is zero-stability. A root outside the unit disk causes geometric amplification as $n\sim T/h$ increases; a repeated unit-circle root permits polynomial growth in $n$. Either can destroy convergence even for starting perturbations tending to zero.

The Dahlquist equivalence theorem states that, for a fixed-coefficient method and a well-posed ordinary differential equation with a suitably Lipschitz continuous right-hand side, consistency of a numerical method and zero-stability are equivalent to convergence for all convergent starting data. If the starting errors are $O(h^p)$ and the exact solution is sufficiently smooth, the global error is $O(h^p)$ on a fixed time interval. A stability estimate bounds the propagated sum of $O(h^{p+1})$ local defects over $O(1/h)$ steps; the discrete Gronwall inequality controls the dependence of $f$ on the accumulated error. Implicit methods also require the locally consistent, uniquely solvable update branch for sufficiently small $h$.

For example, the explicit Euler method has order one and $\rho (\zeta )=\zeta -1$, satisfying the root condition for a multistep method. The two-step Adams-Bashforth method has order two and $\rho (\zeta )=\zeta (\zeta -1)$, again satisfying it. The trapezoidal rule has order two with one step. The two-step Simpson multistep method, $Y_{n+2}-Y_n=h(f_{n+2}+4f_{n+1}+f_n)/3$, has order four and the simple roots $1,-1$: it is zero-stable, although its parasitic amplification root can spoil long-time decay. These examples distinguish finite-time numerical convergence from absolute stability.

The first Dahlquist barrier gives the maximum possible order of a convergent ordinary real $k$-step method:These bounds assume only first-derivative evaluations, fixed coefficients and zero-stability; the sixth-order multiderivative multistep method in solution is outside that class.

To obtain arbitrarily high convergent order within the ordinary class, increase the number of steps. The Adams-Bashforth method integrates a polynomial interpolation of $k$ past derivative values and has order $k$. The Adams–Moulton method includes the new endpoint and has order $k+1$ for $k$ steps. In both families $\rho (\zeta )={\zeta }^{k-1}(\zeta -1)$, so the root condition for a multistep method holds for every fixed $k$. Starting values can be generated by a Runge-Kutta method of sufficient accuracy. The backward differentiation formula instead differentiates an interpolating polynomial; its order-$k$ members are zero-stable only for $1\le k\le 6$.

Finally, high order and favorable stiff stability are separate design goals. By the Second Dahlquist barrier, an ordinary A-stable linear multistep method has order at most two. One may accept a smaller linear stability domain, use suitable backward differentiation formula members for stiff problems, or switch to high-order implicit Runge-Kutta methods or multiderivative multistep methods. A convergent high-order formula still needs time steps appropriate to its stability properties and the problem's eigenvalues.

A spatial finite difference method and time discretization for a linear evolution problem give a recurrence $U^{n+1}=B_{h,k}{U}^n$. To obtain numerical stability on $0\le nk\le T$, one needs a bound $\Vert B_{h,k}^n\Vert \le C_T$ uniform in the permitted grid sizes and time steps. More generally a growth estimate $\Vert B_{h,k}^n\Vert \le C{e}^{cnk}$ is sufficient. The finite-time stability versus power boundedness distinction matters: growth per step of $1+O(k)$ can be compatible with the latter finite-time estimate, whereas unrestricted geometric growth at a fixed Courant number is not.

The eigenvalue method diagonalizes $B_{h,k}$ into modes. If $B_{h,k}$ is a normal matrix, the spectral theorem gives $\Vert B_{h,k}^n{\Vert }_2={\max }_j|{\mu }_j{|}^n$. Thus $|{\mu }_j|\le 1$ gives contraction, and $|{\mu }_j|\le e^{ck}$ gives the preceding growth estimate. A diagonalizable matrix $B=V\Lambda V^{-1}$ has $\Vert B^n{\Vert }_2\le \Vert V{\Vert }_2\Vert V^{-1}{\Vert }_2{\max }_j|{\mu }_j{|}^n$, so its condition number must remain bounded uniformly as the mesh changes. A stable-looking list of eigenvalues alone is insufficient.

Indeed,Both eigenvalues are one, but the Jordan block gives growth proportional to $n\sim T/k$, violating uniform finite-time numerical stability. A nonnormal Fourier amplification matrix can similarly exhibit large amplification even with its eigenvalues inside the unit disk. Eigenvector conditioning and the full powers of the amplification operator are essential.

For constant-coefficient schemes on a periodic or infinite grid, a Fourier transform diagonalizes scalar translation-invariant stencils. This yields von Neumann stability analysis using the Fourier amplification symbol $G(\theta )$. The discrete Parseval identity turns a bound on $G(\theta {)}^n$ into a discrete L2 norm estimate. For systems, the symbol is a matrix; a uniform power bound for matrix Fourier symbols is required, rather than a bound on its spectral radius alone.

For the heat equation, a central finite difference in space and the explicit Euler method in time giveThe scalar condition $|G|\le 1$ for every frequency is exactly $0\le r\le 1/2$. With zero Dirichlet boundary conditions, the spatial matrix is symmetric and diagonalized by sine modes, leading to the same uniform sufficient restriction. By contrast, the Backward Euler diffusion scheme has $G(\theta )=[1+4r{\sin }^2(\theta /2){]}^{-1}$ and is stable for every $r\ge 0$.

For the advection equation $u_t+a{u}_x=0$, $a>0$, centered space with the explicit Euler method gives $G=1-i\nu \sin \theta$, $\nu =ak/h$, and $|G{|}^2=1+{\nu }^2{\sin }^2\theta$. At fixed positive Courant number, repeated powers are unbounded as the mesh is refined. The upwind finite difference scheme instead haswhich gives contraction precisely for $0\le \nu \le 1$. This comparison shows why the spatial stencil and time method must be analyzed together.

Boundary treatment can destroy the normality or mode decomposition available on a periodic grid, so a finite-interval scheme needs its own uniform operator estimate. Multilevel time schemes likewise require the root condition for a multistep method and uniform bounds for the associated amplification matrices; simple scalar root tests do not establish these bounds for arbitrary systems. Once consistency of a numerical method, a well-posed linear evolution problem and the required numerical stability estimate hold, the Lax equivalence theorem gives numerical convergence. Eigenvalues provide a powerful route to that estimate when the accompanying eigenvector and boundary analysis is justified.