Past exam of the mathematics course of the University of Cambridge / 2003 / iii / Paper 67 / 3

There is a genuine indexing defect in the printed PDF: its new value is at $m+1$, while every old-time coefficient is centered at $m$. On a fixed grid the literal residual containsFor a general solution the leading term is $u_x/(\mu \delta )$, so this is not a consistent diffusion discretization. It also leaves one new interior value unspecified while imposing a generally incompatible update on a prescribed boundary value. The calculation that follows is for the intended same-node update $U_m^{n+1}=U_m^n+k{L}_{\delta }{U}_m^n$; that correction is necessary, not an unnoticed transcription change.

For the midpoint flux stencil for one-dimensional diffusion, defineAssume sufficient smoothness of $a$ and the solution for the displayed remainders. Multiplying the centered Taylor expansions of $a$ and $u$ givesThe first time difference is $(u(t+k)-u(t))/k=u_t+k{u}_{tt}/2+O(k^2)$. After using the diffusion equation, the normalized local truncation error isTherefore, for the corrected update,The first is the residual divided by the time step; the second is the unnormalized one-step defect. At fixed $\mu =k/{\delta }^2$ they are respectively $O({\delta }^2)$ and $O({\delta }^4)$. A rate based on derivatives is not justified for a merely bounded nonsmooth $a$.

Use the corrected same-node method identified above, with $\delta =1/M$, interior indices $1,\dots ,M-1$ and zero endpoint values. Its spatial matrix $L_{\delta }$ is real symmetric. Summation by parts gives, for a vector extended by zero at the endpoints,Hence $-L_{\delta }$ is positive definite, and ${\sum }_j(v_{j+1}-v_j{)}^2\le 4{\sum }_j{v}_j^2$ bounds its eigenvalues by $4a_{+}/{\delta }^2$. The forward Euler method has amplification matrix $I+k{L}_{\delta }$. Its eigenvalues lie in $[1-4\mu a_{+},1]$, so their moduli are at most one whenSymmetry makes this a contraction in the mesh-weighted discrete L2 norm, uniformly over all grids and all coefficient samples satisfying the bounds. There is also a maximum norm proof: the update weights are $\mu a_{m-1/2}$, $1-\mu (a_{m-1/2}+a_{m+1/2})$, and $\mu a_{m+1/2}$. They are nonnegative and sum to one, giving a contraction after the boundary values are included.

This mesh-independent bound is sharp. For $a\equiv a_{+}$, the spatial eigenvectors are $\sin (j\pi m/M)$ and their amplification factors areIf $\mu >1/(2a_{+})$, sufficiently fine grids have a highest-frequency factor below $-1$, giving exponential growth in the step number.

On one fixed grid, rather than uniformly under refinement, the sharp universal stability of a numerical method bound in the discrete L2 norm is slightly larger:Indeed the quadratic form above is bounded by $a_{+}$ times the constant-coefficient difference form, whose largest eigenvalue is $4{\cos }^2(\pi /(2M))/{\delta }^2$; equality is realized by the constant coefficient. Its limit as $M\to \infty$ is the boxed bound. No such fixed-grid matrix conclusion repairs the ill-defined shifted-boundary update in the literal printed formula.