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

Write the Runge-Kutta method coefficients as $A\in {\mathbb{R}}^{\nu \times \nu }$, $b\in {\mathbb{R}}^{\nu }$, and let $e$ be the all-ones column. On the test equation $y^{\prime }=\lambda y$, with $z=h\lambda$, the stage vector satisfies $Y=e{y}_n+zAY$ and the next value is $y_{n+1}=y_n+z{b}^{T}Y$. Hence its stability function isBoth polynomials have degree at most $\nu$, and the denominator has value one at zero. Order of a Runge-Kutta method $2\nu$ implies that this test-equation approximation agrees with $e^z$ through degree $2\nu$:By uniqueness of the diagonal Padé approximant, $R=[\nu /\nu {]}_{e^z}$. Explicitly,The permitted diagonal Padé approximant property says this rational function has no poles in the closed left half-plane and satisfies $|R(z)|\le 1$ there. ThereforeThere is also no hidden stage-solve singularity from a canceled denominator: the Padé approximant numerator and denominator are coprime and each has degree $\nu$. The original determinant denominator already has degree at most $\nu$, so, with its normalization at zero, it equals $P_{\nu }(-z)$ and cannot contain an additional canceled factor. This proves maximal-order Runge-Kutta methods are A-stable.

Use the Gauss collocation coefficient construction. Find the $\nu$ polynomial roots $c_1,\dots ,c_{\nu }$ of the shifted Legendre polynomial $P_{\nu }(2c-1)$ in $(0,1)$, and formThen setThese are explicit polynomial integrals. For example, if ${\ell }_j(s)={\sum }_{r=0}^{\nu -1}{d}_{jr}{s}^r$, calculate $a_{ij}={\sum }_r{d}_{jr}{c}_i^{r+1}/(r+1)$ and $b_j={\sum }_r{d}_{jr}/(r+1)$. The resulting implicit Gauss collocation method has order $2\nu$; its Gaussian quadrature nodes integrate polynomials through degree $2\nu -1$ exactly. The one-stage example is the implicit midpoint rule, with $c_1=a_{11}=1/2$ and $b_1=1$. No search over nonlinear Runge-Kutta order conditions systems is needed.

Put $h=\Delta x$, $k=\Delta t$, with fixed positive $\mu =k/h$. Insert a smooth exact solution and divide the residual by $2k$:Taylor expansion givesThe differential equation cancels the leading part, leaving normalized local truncation error $O(k^2+h^2)$. Since $u_{ttt}=({\partial }_x+{\partial }_y{)}^{3}u$ contains mixed derivatives, there is no fixed positive Courant number canceling this leading error for every smooth solution. Thus the method is second order in time and space. Global second-order convergence also requires stability and initial values at both time levels accurate enough to retain this order.

Use a spatial discrete Fourier mode $u_{k,j}^n=a_n{e}^{i(k\xi +j\eta )}$. Its recurrence and amplification polynomial of a multilevel finite difference scheme areThe polynomial roots are $G_{\pm }=is\pm \sqrt{1-s^2}$. For $|s|<1$ both have modulus one and are separated. Since $|\sin \xi +\sin \eta |\le 2$, a fixed $0<\mu <1/2$ bounds their separation below by $2\sqrt{1-4{\mu }^2}$. The two-level companion matrix is therefore diagonalizable with uniformly bounded eigenvector matrix and inverse. To see the uniform bound explicitly, write $V=\left(\begin{array}{cc}{G}_{+}& G_{-}\\ 1& 1\end{array}\right)$. Its determinant is $G_{+}-G_{-}$; all its entries have modulus one, and those of $V^{-1}$ are bounded by the reciprocal of the root separation. The companion matrix powers are $V\mathrm{diag}(G_{+}^n,G_{-}^n)V^{-1}$, bounded independently of frequency and step number. Parseval's identity converts this into a mesh-independent discrete $L^2$ bound for arbitrary perturbations at the two initial levels.

If $\mu >1/2$, choose $\xi =\eta =\pi /2$. Then $s=2\mu >1$, and one polynomial root has modulus $s+\sqrt{s^2-1}>1$. On periodic meshes with the number of points divisible by four, this is an actual grid mode, producing exponential linear instability.

At $\mu =1/2$ the same phases give $(G-i{)}^2$. The companion matrix is not a scalar matrix, so the double polynomial root has a nontrivial Jordan block. The solution $a_n=n{i}^n$ has bounded starting amplitudes but grows like $n$. On a fixed physical time interval $n$ is of order $1/k$, so no mesh-independent stability bound exists. Consequentlyis the stability interval under the standard two-level definition. Testing only polynomial root moduli would misleadingly include the endpoint. A special startup selecting the non-growing branch at that endpoint restricts the perturbations and does not establish stability of the full recurrence. This is the two-dimensional leapfrog stability threshold.

The first characteristic polynomial factors asAlso $\rho (1)=0$ and ${\rho }^{\prime }(1)=\sigma (1)=2(1-\alpha )$, so the order conditions for a linear multistep method hold for every $\alpha$. Convergence of the full linear multistep method recurrence requires, in addition, the root condition for a multistep method: all polynomial roots of $\rho$ lie in the closed unit disk and its polynomial roots on the unit circle are simple.

For $-1<\alpha <1$, write $\alpha =\cos \theta$ with $0<\theta <\pi$. The other polynomial roots are $e^{\pm i\theta }$, distinct from each other and from 1, so the root condition for a multistep method holds. For $|\alpha |>1$, the two quadratic polynomial roots are real and reciprocal, and one has modulus greater than one. At $\alpha =1$, $\rho =(w-1{)}^3$; at $\alpha =-1$, $\rho =(w-1)(w+1{)}^2$. Both endpoints violate simplicity. Consistency of a numerical method plus zero-stability therefore givesThe conclusion applies to the recurrence as printed, including its starting data. Canceling common factors is not innocuous: at $\alpha =1$ the polynomials contain $(w-1{)}^2$, and at $\alpha =-13/5$ they contain $w+5$. The reduced recurrences exclude parasitic solutions of the original one. This is the common-factor cancellation defect in a multistep recurrence, so those parameter values are not added to the interval.

Apply the exponential-symbol order criterion for a multistep method, which follows by substituting the bilateral shift operator $e^{h{\partial }_t}$ into the exact-solution residual. Direct expansion givesAll coefficients through degree three vanish. If $\alpha \ne -5$, the fourth-degree coefficient is nonzero, so the normalized local truncation error has exact order three. At $\alpha =-5$, the fourth-degree coefficient vanishes and the fifth-degree coefficient is $1/10$, which is nonzero. ThusThis is the formal order of the full recurrence. The order-four member is not zero-stable, so it is not convergent. All the convergent members from part (a) have order three. At the degenerate value $\alpha =1$, the full recurrence still has formal order three, while cancellation of its common factor gives Backward Euler method of order one; that cancellation changes the method rather than its order calculation.

For $|\alpha |>1$, an exterior polynomial root of $\rho$ persists by continuity for sufficiently small negative real $z$ in the stability polynomial $\rho (w)-z\sigma (w)$. Hence these parameters cannot be A-stable.

For $-1\le \alpha <1$, the principal polynomial root at $w=1,z=0$ is simple. Let $s(z)=\log w(z)$ denote its analytic logarithm near zero. Using the expansion in part (b), and ${\rho }^{\prime }(1)=2(1-\alpha )$, solve the polynomial root equation to obtainIndeed the residual at $s=z$ is $-(\alpha +5)z^4/12$, and the correction in $s$ cancels it by multiplication with $2(1-\alpha )$. Choose $z=-C{\epsilon }^4/2+i\epsilon$. This lies strictly in the left half-plane, butfor small positive $\epsilon$. Thus $|w(z)|=e^{\mathrm{Re}s(z)}>1$, proving linear instability within the A-stability domain. This is a direct principal-root obstruction to A-stability, rather than an appeal to an order-barrier theorem alone.

Finally, at $\alpha =1$, the full stability polynomial isIt retains a double unit polynomial root for every negative $z$, allowing growing parasitic solutions. ThereforeThe canceled Backward Euler method recurrence at $\alpha =1$ is A-stable, but its admissible solutions and starting relations differ from those of the original method.

Take $L=-d^2/d{x}^2+x$. The Airy equation is $Lu=0$. To discuss positivity of this linear operator, use its homogeneous variation domain $D(L)=\{v\in H^2(0,1):v(0)=0,v^{\prime }(1)=0\}$. Integration by parts givesThe boundary term is zero at both ends. The same calculation with two different functions proves that the associated form is a symmetric bilinear form. Positivity here concerns the homogeneous domain; the condition $u(0)=1$ instead defines an affine set of admissible solutions.

For the variational problem defineand minimizeThe Sobolev trace at zero is meaningful in $H^1$. No derivative boundary value is imposed on this trial space: the Neumann boundary condition will emerge naturally. For $v\in V$, Cauchy-Schwarz inequality and integration giveso $a(v,v)\ge \frac{2}{3}\Vert v{\Vert }_{H^1}^2$. Thus $a$ is a bounded, symmetric, coercive bilinear form on the Hilbert space $V$. The Lax-Milgram theorem says that a bounded coercive bilinear form and a bounded linear functional determine a unique weak solution. Apply it toand set $u=1+v$. This is precisely the first variation condition $a(u,w)=0$ for $J$.

Compactly supported test functions first yield $u^{\prime \prime }=xu$ in the sense of distributional derivatives. Since $xu\in L^2$, the solution is in $H^2$. Integration by parts then leaves $u^{\prime }(1)w(1)=0$ for every $w\in V$, hence $u^{\prime }(1)=0$. The other boundary condition is built into $1+V$. Conversely a solution of this boundary value problem satisfies the first variation condition. Finally,This proves existence, uniqueness and the global minimum characterization, including the natural boundary condition. It is the mixed-boundary Airy energy principle.

Choose a partition $0=x_0<x_1<\cdots <x_N=1$, with the usual nodal piecewise-linear hat functions ${\varphi }_i$, including the half hats at the endpoints. Write $u_h={\sum }_{i=0}^N{U}_i{\varphi }_i$ with $U_0=1$. Differentiating the Ritz method energy with respect to $U_i$, $1\le i\le N$, givesOn an element $[a,b]$ of length $\ell$, its two basis functions are $(b-x)/\ell$ and $(x-a)/\ell$. Direct integration yields the element matrixFor example the first weighted diagonal is ${\int }_a^{b}x(b-x{)}^2/{\ell }^{2}dx=\ell (3a+b)/12$; the off-diagonal is ${\int }_a^{b}x(b-x)(x-a)/{\ell }^{2}dx=\ell (a+b)/12$. This gives the affine-weighted hat mass matrix and explicit equations on any partition by assembling adjacent elements.

For the uniform choice $h=1/N$, $x_i=ih$, the assembled interior equations areThe last node contributes only one element, so its equation isIn the first interior equation the known $U_0$ term moves to the right, giving $1/h-h^2/12$. For $N=1$ only the endpoint equation is needed. The last row enforces the weak Neumann boundary condition, so no artificial value outside the interval is introduced. The matrix on $U_1,\dots ,U_N$ is a symmetric positive-definite matrix: its quadratic form is $a(v_h,v_h)>0$ for every nonzero trial variation $v_h$. The discrete minimizer is therefore unique.

Let $h=\Delta x$, $k=\mu h$ and $M_r={\sum }_{j\in \{-2,-1,1,2\}}{a}_j{j}^r$. Along a smooth solution, ${\partial }_t^{r}u={\partial }_x^{r}u$. The exact difference $u(x,t+k)-u(x,t-k)$ contains only odd derivatives. Matching terms through degree four therefore requiresThe first two conditions give $a_{-1}+a_1=a_{-2}+a_2=0$, which also implies $M_4=0$. With $a_{-j}=-a_j$, the remaining equations are $a_1+2a_2=\mu$ and $a_1+8a_2={\mu }^3$. HenceThe fifth spatial moment is $M_5=10{\mu }^3-8\mu$. After dividing the exact-minus-scheme residual by $2k$, its first possible nonzero term isThis is the fourth-order two-step advection stencil. Thus the normalized local truncation error is $O(h^4)$ at fixed positive Courant number, as required. At $\mu =1$ or $2$ the formula translates the physical advection branch exactly, but consistency of a numerical method alone does not decide stability of all two-level perturbations.

Fourier transform on the infinite spatial lattice giveswith polynomial roots $G_{\pm }=iA\pm \sqrt{1-A^2}$. When $|A|$ is uniformly less than one, the polynomial roots have modulus one and a uniform positive separation. Diagonalizing the companion matrix then bounds all its powers independently of $\theta$ and $n$. Parseval's identity proves stability in the discrete $L^2$ norm for arbitrary square-summable starting perturbations at both levels.

At $\mu =1/2$ the coefficients giveConsequently $\mu =1/2$ is stable.

At $\mu =3/2$,The growing polynomial root at this phase has modulusTo turn this into a Cauchy problem linear instability proof, rather than rely on a plane wave of infinite norm, choose a small interval around $\pi /3$ on which the growing polynomial root has modulus at least $1+\delta$. Choose an initial Fourier transform amplitude that is a nonzero square-integrable function supported there, and set the level-one amplitude equal to that polynomial root times the level-zero amplitude. The recurrence then multiplies by the same growing polynomial root at every step on the support. Parseval identity gives $\Vert u^n{\Vert }_2\ge (1+\delta {)}^n\Vert u^0{\Vert }_2$, whereas both initial norms are bounded independently of $h$. Real data can be obtained by adding the conjugate interval around $-\pi /3$. Since $n=T/(\mu h)$ tends to infinity at fixed $T$, $\mu =3/2$ is unstable.

The Engquist-Osher method is an upwind finite volume method. For cell values $U_j^n$, mesh width $h$ and time step $k$, put $r=k/h$ and use the conservative updatewhere the numerical flux isConsistency of a numerical method follows from $F(u,u)=f(u)$. Positive characteristic speeds use the left state, and negative speeds use the right state. The flux difference telescopes when summed over cells, giving a discrete conservation law. The method is first order in space and time for smooth solutions; its conservation law form is also appropriate when shock waves develop.

Let $s_{\ast }$ be the unique sonic point of a scalar flux, $f^{\prime }(s_{\ast })=0$. Convexity implies $f^{\prime }<0$ to its left and $f^{\prime }>0$ to its right. Put $f_{\ast }=f(s_{\ast })$ and splitThen $f=f_{\ast }+f^{+}+f^{-}$, $(f^{+}{)}^{\prime }\ge 0$, $(f^{-}{)}^{\prime }\le 0$, andThis is the sonic-point splitting of a convex Engquist-Osher flux. If both states are above $s_{\ast }$, $F=f(a)$; if both are below it, $F=f(b)$. Across a transsonic rarefaction wave, $a\le s_{\ast }\le b$, the flux is the minimum $f_{\ast }$. In the opposite crossing it is $f(a)+f(b)-f_{\ast }$, which in general differs from the exact Godunov numerical flux. This distinction does not affect the monotonicity proof.

Here stability requires the usual Courant–Friedrichs–Lewy condition. Let the initial values lie in $[m,M]$, and chooseA differentiable convex function has a continuous derivative on this interval, so the maximum is finite. Write the three-input update as $H(a,b,c)$. Differentiating givesThus the update $T$ preserves componentwise order. Constants are fixed, since $H(c,c,c)=c$, so comparison with the constants $m,M$ proves $m\le U_j^n\le M$ for every time level. In particular the same CFL condition remains applicable. This is a mesh-independent maximum-principle stability bound.

A stronger stability statement is contraction in discrete $L^1$. Consider two solutions $U,V$ in the same invariant interval, on a periodic mesh or on the infinite lattice with summable differences. Let $Q=U\vee V$ be their componentwise maximum. Monotonicity gives $TQ\ge TU,TV$, henceThe discrete conservation law for the difference gives the equality: on a finite periodic mesh the flux sum cancels exactly, and on the infinite lattice it follows by truncation and the fact that the flux is a Lipschitz function on $[m,M]$. Interchanging $U,V$ and adding provesThis is L1 contraction of a monotone conservative scheme. Since $T$ commutes with the cell shift, applying this contraction to $U$ and its shift gives ${\sum }_j|U_{j+1}^{n+1}-U_j^{n+1}|\le {\sum }_j|U_{j+1}^n-U_j^n|$. The scheme is therefore a total variation diminishing scheme as well.

The same order argument proves a discrete entropy inequality, explaining why this stable discretization selects entropy solutions rather than nonphysical expansion shock waves. For a constant $c\in [m,M]$, setBecause $T(U\vee c)\ge TU,c$ and $T(U\wedge c)\le TU,c$, their componentwise difference is at least $|TU-c|$. Subtracting their updates givesFor constants outside $[m,M]$, the entropy inequality is an equality by the invariant bound and the conservative update. At equal states $Q_c(u,u)=\mathrm{sgn}(u-c)[f(u)-f(c)]$, the Kruzhkov entropy flux. The flux is a Lipschitz function on the invariant interval: writing $L=\max |f^{\prime }|$, the conservative update gives $h{\sum }_j|U_j^{n+1}-U_j^n|\le 2kL{\sum }_j|U_{j+1}^n-U_j^n|$. This time-translation bound, the invariant bound and the total variation estimate give local space-time compactness for bounded-variation initial conditions as the mesh is refined; the conservative update passes to the weak conservation law, and the displayed inequality passes to its Kruzhkov entropy inequality. This identifies the limit as the entropy solution.

Finally the time-step qualification is essential. For the convex Inviscid Burgers equation flux $f(u)=u^2/2$, whose unique sonic point of a scalar flux is zero, linearize about a positive constant $U$. The perturbation update becomes $v_j^{n+1}=(1-rU)v_j^n+rU{v}_{j-1}^n$. At Fourier frequency $\pi$ its multiplier is $1-2rU$, with modulus greater than one if $rU>1$. Wave packets near this phase give linear instability even for summable Cauchy perturbations. Thus the assumptions on $f$ imply the stability results above under the CFL condition, not for arbitrary $k/h$. This is CFL necessity for explicit Engquist-Osher stability.

A Mehrstellen method gains accuracy by exploiting the differential equation in the local truncation error of a compact finite difference method stencil. The extra nearby values alone do not guarantee higher order: one must also modify the source term. Consider $\Delta u=f$ on a square with prescribed Dirichlet boundary data, using the same spacing $h$ in both directions.

For comparison, the five-point Laplacian has expansionIts leading error is not a multiple of ${\Delta }^{2}u$, because the mixed fourth derivative is missing. To obtain an isotropic leading error, consider a symmetric nine-point finite-difference stencil with weight $a$ on each axial neighbor, $b$ on each diagonal neighbor, and $c$ on the center, all divided by $h^2$. Vanishing on constants and consistency of a numerical method require $4a+4b+c=0$ and $a+2b=1$. Taylor expansion then gives the fourth-derivative contributionFor this to equal $h^2{\Delta }^{2}u/12$ we must choose $b=1/6$, hence $a=2/3$, $c=-10/3$. The resulting nine-point finite-difference stencil operator isIts expansion, retaining the next term for later use, isThus $D_9$ by itself is second order on a general function. On a Poisson equation solution, however, ${\Delta }^{2}u=\Delta f$. Approximating this known source derivative by $D_{5}f$ gives the fourth-order Mehrstellen method equationEquivalently, its two stencils areThe normalized local truncation error is $O(h^4)$, since $D_{5}f-\Delta f=O(h^2)$. If the equation is multiplied by $6h^2$, the row residual is $O(h^6)$; this change of normalization does not make the solution sixth order. The fourth-order correction of the nine-point Poisson stencil retains only the closest axial and diagonal unknowns. A fourth-order one-dimensional second derivative instead uses $(-U_{i+2}+16U_{i+1}-30U_i+16U_{i-1}-U_{i-2})/(12h^2)$, widening the unknown stencil and introducing negative outer weights. Mehrstellen method accuracy is obtained by correcting known source values rather than widening that unknown stencil.

The compact scheme has useful stability and solvability properties. With homogeneous Dirichlet boundary data, $-D_9$ is symmetric, has positive diagonal and nonpositive off-diagonal entries, and is a positive definite symmetric operator. Indeed its quadratic form is a sum of positive edge weights times squared differences, with boundary values set to zero. Vanishing of this sum would force a constant on the connected grid and then zero by connection to the boundary. Consequently the Dirichlet boundary condition linear system has a unique solution. The same weights give a discrete maximum principle: if $D_{9}v\ge 0$ at every interior node and $v\le 0$ on the boundary, a positive maximum inside would make every weighted difference to its neighbors nonpositive. Equality forces the maximum to propagate to the boundary, which is impossible. Thus $v\le 0$ everywhere.

This also turns consistency of a numerical method into a quantitative convergence result, rather than merely suggesting an order. Let $B_{h}f=f+h^2{D}_{5}f/12$, let $\tau =D_{9}u-B_{h}f$ be the normalized residual of the exact solution, and set $e=U-u$. Exact Dirichlet boundary data imply $e=0$ on the boundary and $D_{9}e=-\tau$. On the unit square the quadratic barrierhas $D_{9}q=-1$ exactly, is nonnegative on the boundary, and has maximum $1/8$. If $\Vert \tau {\Vert }_{\infty }\le \epsilon$, then $D_9(e-\epsilon q)\ge 0$ and $D_9(-e-\epsilon q)\ge 0$, while both comparison functions are nonpositive on the boundary. The discrete maximum principle provesThis is maximum-norm convergence of the corrected nine-point Poisson scheme. It requires sufficient smoothness of the solution up to the boundary for the uniform truncation bound. On other domains the source and boundary closures must retain the required accuracy; boundary singularities can invalidate a smooth-solution order estimate.

The same expansion reveals how to add still more accuracy. Differentiating $\Delta u=f$ givesThus the sixth derivatives appearing in the $h^4$ term can also be expressed solely through the source. With exact source derivatives,is the sixth-order source correction of the nine-point Poisson stencil, with normalized residual $O(h^6)$. If source derivatives are replaced by differences, $\Delta f$ must be approximated through fourth order and the fourth derivatives through second order to retain that accuracy. The unknowns still occupy the same nine-point finite-difference stencil; only evaluation of known source data becomes more elaborate. The same inverse bound yields sixth-order convergence when these approximations and the boundary data have matching accuracy.

A particularly instructive case is the Laplace equation, $f=0$. Both displayed error corrections automatically vanish: ${\Delta }^{2}u=0$ and $u_{xxxxyy}+u_{xxyyyy}=\Delta u_{xxyy}=0$, which also makes the pure sixth derivatives sum to zero. Expanding the stencil two orders further gives, for a sufficiently smooth harmonic function,For verification, the three eighth-derivative contributions are $(u_{xxxxxxxx}+u_{yyyyyyyy})/20160$, $(u_{xxxxxxyy}+u_{xxyyyyyy})/2160$ and $u_{xxxxyyyy}/864$. Because $u$ is a harmonic function, their sum is $u_{xxxxxxxx}(1/10080-1/1080+1/864)=u_{xxxxxxxx}/3024$. The example $u=\mathrm{Re}(x+iy{)}^8$ has $D_{9}u(0,0)=40h^6/3$, so this term is genuinely present. The unmodified harmonic nine-point scheme therefore has sixth-order nodal convergence under the same smoothness and boundary assumptions. This harmonic superconvergence of the nine-point stencil is special to harmonic functions; it should not be confused with the fourth-order source-corrected scheme for general Poisson equation data.

The appeal of Mehrstellenverfahren is this combination of a compact sparse matrix, favorable discrete maximum principle properties, and high accuracy obtained by using the governing equation. The derivation also shows precisely which accuracy belongs to the operator on arbitrary functions, which belongs to its action on solutions, and which survives in the computed boundary-value solution.