Numerical solutions of partial differential equations course of the University of Oxford (B6.1)

Original solutions to past examinations of the numerical solutions of partial differential equations course of the University of Oxford, using finite difference methods, energy estimates and Fourier stability analysis.

Original solutions to all three questions of the numerical solutions of partial differential equations course of the University of Oxford Part B examination, held on 6 June 2023. Source: original released PDF, bundle pages 292–295; paper pages 1–4. The numbering follows the paper. In 3(e), the growth constant must be allowed to depend on the mesh; the literal mesh-independent assertion in the native PDF is false. The corrected result and the intended failure of practical stability of a finite difference scheme are proved below.

Write $w=u-v$ and $r=f-g$. Subtraction of the two boundary value problems gives $-w^{\prime \prime }+w+u^3-v^3=r$, with $w(0)=w(1)=0$. Multiply by $w$ and apply integration by parts:Consequently the required energy identity isThe nonlinear term is nonnegative pointwise, sinceThis is the monotonicity energy estimate for a semilinear Dirichlet problem for the monotone function $s\mapsto s^3$. Apply Young's inequality for products, $rw\le (r^2+w^2)/2$, inside the last integral and discard the nonnegative nonlinear term. Moving $\frac{1}{2}{\int }_0^1{w}^{2}dx$ to the left givesIn particular, equal forcing gives $w=0$, proving uniqueness and a quantitative continuous dependence on forcing estimate for this boundary value problem.

Let $f_i=f(ih)$ and $g_i=g(ih)$, and define the scaled forward difference operator, backward difference operator and second-order central difference byThe two finite difference methods are the nonlinear systemswith the prescribed Dirichlet boundary conditions $U_0=V_0=A$ and $U_N=V_N=B$. These equations determine the interior nodal values.

For clarity distinguish the interior and edge inner products and their norms:The bracketed norm includes the last edge $i=N$. Let $W=U-V$, so $W_0=W_N=0$, and put $q_i=D_x^{-}{W}_i$ for $1\le i\le N$. Direct telescoping gives discrete summation by parts:Subtract the two nonlinear systems, multiply each equation by $h{W}_i$ and sum. The resulting discrete energy identity isThe nonlinear summands are nonnegative by the factorization in (a). Using Young's inequality for products at every node givesThis also proves uniqueness of each discrete solution. Existence can be seen by minimizing, with the endpoints fixed, the strictly convex energy functionalIt is a coercive function of its finitely many interior values, so it attains a minimum. Its first derivatives are $h(-D_x^{+}{D}_x^{-}{W}_i+W_i+W_i^3-f_i)$; thus its minimizer solves the scheme. The same argument applies to $g$.

Use $u_i=u(x_i)$ for the sampled exact solution, and let $e_i=U_i-u_i$. The Dirichlet boundary conditions imply $e_0=e_N=0$. The definition of the truncation error $\phi$ gives the exact discrete equationSubtracting it from the numerical equation therefore givesThe sign on the right is positive because ${\phi }_i$ was defined as the difference approximation minus the exact second derivative. Multiply by $h{e}_i$, sum over the interior nodes and apply discrete summation by parts:Again $(U_i^3-u_i^3)e_i=e_i^2(U_i^2+U_i{u}_i+u_i^2)\ge 0$. By Young's inequality for products, $(\phi ,e{)}_h\le \frac{1}{2}\Vert \phi {\Vert }_h^2+\frac{1}{2}\Vert e{\Vert }_h^2$. HenceNo fourth derivative is needed for this step: the assumed twice continuous differentiability suffices to define $\phi$ and obtain the discrete energy estimate.

On the vector space of nodal values with zero endpoints define the discrete Sobolev normIt is a norm because it is the Euclidean norm of the scaled list of nodal values and first finite differences, and vanishes only for the zero list. It is the discrete counterpart of the Sobolev space $H^1$ norm. Here $W=U-u$ has zero endpoints, so this definition applies directly.

Set $M_4={\max }_{0\le x\le 1}|u^{(4)}(x)|$, finite by the assumed continuity. Taylor's theorem with fourth-order remainders gives, for suitable ${\xi }_{i,+}\in (x_i,x_i+h)$ and ${\xi }_{i,-}\in (x_i-h,x_i)$,Adding the two expansions cancels the odd derivatives, soConsequently $\Vert \phi {\Vert }_h^2\le h(N-1)M_4^2{h}^4/144\le M_4^2{h}^4/144$. Doubling the bound in (c) gives $2\Vert D_x^{-}e]{\Vert }_h^2+\Vert e{\Vert }_h^2\le \Vert \phi {\Vert }_h^2$, and thereforeThus $M_4/12$ is a valid nonnegative constant. To satisfy the request for a strictly positive constant even when $u^{(4)}\equiv 0$, take $C=\max \{1,M_4\}/12$, independent of $h$. This proves second-order convergence of a numerical method in a discrete Sobolev norm, including control of the backward difference operator applied to the error.

For $u^{(m,n)}(x,y)=\sin (m\pi x)\sin (n\pi y)$, direct differentiation yieldsSince the planar Laplacian is $\Delta ={\partial }_{xx}+{\partial }_{yy}$, the required eigenvalue of $-\Delta +I$ isOn the sides $x=0,1$ the first sine factor vanishes, and on $y=0,1$ the second sine factor vanishes. Thus the homogeneous Dirichlet boundary condition holds on the whole boundary. The eigenfunction is nonzero, for instance at $(x,y)=(1/(2m),1/(2n))$, where its value is one.

The five-point Laplacian on the square grid isIt adds the second-order central differences in the two coordinate directions and uses the given zero Dirichlet boundary conditions when a neighbouring node lies on the boundary.

The sine addition formula impliesApply this identity with $\alpha =m\pi h$ in the $x$ direction and with $\alpha =n\pi h$ in the $y$ direction. For $U_{i,j}^{(m,n)}=\sin (m\pi ih)\sin (n\pi jh)$, it givesHence the discrete eigenvalues areBoth sine factors vanish on their respective boundary sides. Also $U_{1,1}^{(m,n)}>0$, since $m\pi /N$ and $n\pi /N$ lie strictly between zero and $\pi$. Thus each array is a nonzero discrete eigenvector, rather than merely the zero solution of the homogeneous system. The separated structure is also described by the five-point Dirichlet Laplacian as a Kronecker sum.

Set $z=\pi h/2$, so $0<z\le \pi /4$. The formulas in (a) and (b) giveFor $0\le t\le \pi /4$, $0\le \sin t\le t$ and $1-\cos t=2{\sin }^2(t/2)\le t^2/2$. ThereforeThese trigonometric inequalities implySubstitution yields the explicit eigenvalue error boundThis constant is asymptotically sharp: Taylor's theorem gives ${\sin }^{2}z=z^2-z^4/3+O(z^6)$, so $({\lambda }_{1,1}-{\Lambda }_{1,1})/h^2\to {\pi }^4/6$. Here $O$ denotes Big O notation.

For every admissible $m,n$, the trigonometric inequality $\sin z\le z$ for $z\ge 0$ givesThus ${\Lambda }_{m,n}\le {\lambda }_{m,n}$, in particular for $m=n=N-1$ and every $N\ge 3$.

For these highest-index eigenvalues, use $h=1/N$ and $\sin ((N-1)\pi /(2N))=\cos (\pi /(2N))$ to obtainDivide by $N^2$ and take the limit:For example, eventually the normalized difference exceeds $(2{\pi }^2-8)/2$, so the unnormalized difference is at least $({\pi }^2-4)N^2$ and tends to $+\infty$. This does not contradict (c): there the indices stay fixed as the grid is refined; here the eigenfunctions oscillate increasingly quickly, with wavelength comparable to the mesh spacing. It illustrates the high-frequency eigenvalue error of the discrete Laplacian.

Use the Backward Euler method in time and the prescribed central finite differences in space. With $j\in \mathbb{Z}$ and $m\ge 0$, the implicit finite difference method for the advection-diffusion equation isEquivalently, for $\nu =\Delta t/\Delta x$ and $\mu =\Delta t/(\Delta x{)}^2$, the new time level solvesThe spatial terms are evaluated at the new time level, making this an implicit scheme. On the whole line there are no finite-endpoint boundary conditions; in the next part its solution is chosen in the l2 sequence space.

For a grid sequence define its semidiscrete Fourier transform in the dimensionless frequency variable $\theta \in [-\pi ,\pi ]$ byFor general square-summable data the transform is interpreted as an limit in the L2 norm, rather than an absolutely convergent pointwise sum. The Parseval identity in this normalization isThe Fourier symbols of $D_x^0$ and $D_x^{+}{D}_x^{-}$ are $i\sin \theta /\Delta x$ and $-4{\sin }^2(\theta /2)/(\Delta x{)}^2$. Transforming the implicit equation givesIts amplification factor is consequentlyThe denominator has real part at least one because $a\ge 0$. It never vanishes, so multiplication by $G_I$ defines a unique next value in the l2 sequence space, andThe Parseval identity now givesThis is unconditional practical stability of a finite difference scheme: in the course's terminology the numerical evolution does not increase this norm. Linearity gives the same estimate for perturbations of the initial data. This contraction convention is stronger than a finite-time mesh-uniform bound allowing physical exponential growth.^[ref]

Use the forward Euler method in time and evaluate both spatial finite difference operators at the old time level:With $\nu =\Delta t/\Delta x$ and $\mu =\Delta t/(\Delta x{)}^2$, the explicit update isEvery value on the new time level is computed directly from three old values. The sum of the stencil weights is one, as required to preserve constant solutions of the advection-diffusion equation.

Using the semidiscrete Fourier transform from (b), the explicit scheme has amplification factorSet $r=a\mu$ and $s={\sin }^2(\theta /2)\in [0,1]$, so ${\sin }^2\theta =4s(1-s)$. Direct expansion givesThe expression in square brackets is an affine function of $s$. At the two endpoints its values areusing precisely ${\nu }^2\le 2r\le 1$. Every intermediate value is a convex combination of these endpoints, so $|G_E{|}^2\le 1$ for every frequency. The Parseval identity provesThus the scheme is practically stable under the given condition. For $a>0$ it is equivalent to $\Delta t\le 2a$ and $\Delta t\le (\Delta x{)}^2/(2a)$. This is the Forward Euler stability for centered advection-diffusion calculation with the diffusion coefficient retained. The proof does not require all three stencil weights to be nonnegative.

When $a=0$, the explicit amplification factor becomesA single Fourier mode has constant modulus in space and is not in the l2 sequence space. Instead choose a whole interval of frequencies on which there is definite growth. LetThe inverse Fourier transform here gives the real sequenceBy the Parseval identity, its squared norm is $\Vert U^0{\Vert }_{{\ell }_2}^2=\Delta x|E|/(2\pi )=\Delta x/2$, finite and strictly positive. These data are also compatible with sampling a bounded continuous function: useThe value at zero is the limit of the displayed quotient, and $u_0(j\Delta x)=U_j^0$.

On $E$ we have ${\sin }^2\theta \ge 1/2$. Therefore the transformed solution and the Parseval identity giveFor each fixed pair $\Delta x,\Delta t>0$, takeTaking square roots proves $\Vert U^m{\Vert }_{{\ell }_2}\ge (1+c{)}^m\Vert U^0{\Vert }_{{\ell }_2}$ for every $m\ge 1$. In particular the norm already increases at the first step, so no positive mesh pair gives practical stability in the contraction sense of (b).

The native PDF places the existence of $c>0$ before the universal quantifier over mesh pairs. If this means a single positive $c$ independent of $\Delta x,\Delta t$, that assertion is false for nonzero data. Indeed, for every datum the same Fourier symbol givesFor any proposed fixed $c>0$, choose $0<\nu <\sqrt{(1+c{)}^2-1}$. The upper bound then contradicts the requested lower bound at $m=1$. The correct quantifiers are $\forall \Delta x,\Delta t>0\exists c=c(\Delta x,\Delta t)>0$, with nonzero square-summable data as constructed above.

Finally, failure of this stronger practical stability of a finite difference scheme should not be confused with every possible refinement notion of numerical stability. If $\Delta t/(\Delta x{)}^2\le K$, then $m\Delta t\le T$ givesThus a finite-time mesh-uniform bound can hold under this severe refinement even though contraction fails on every positive mesh. This is the distinction recorded in centered-advection refinement-dependent stability.