Past exam of the mathematics course of the University of Cambridge / 2008 / iii / Paper 80 (Numerical Solution of Partial Differential Equations)

Apply the exponential-symbol order criterion for a multistep method. The exact-solution residual of the linear multistep method is encoded by $E(z)=\rho (e^z)-z\sigma (e^z)$; vanishing through degree $p$ is equivalent to the order conditions for a linear multistep method through order $p$. Here $\rho (1)=0$, ${\rho }^{\prime }(1)=\sigma (1)=1-\alpha$, and direct Taylor expansion givesFor example, for $n\ge 2$ the coefficient of $z^n$ can be calculated asFor $n=2,3$ this is zero, verifying the lower conditions explicitly. The constant and linear coefficients also vanish for every $\alpha$. Therefore the formal order is at least three for every $\alpha$.

The fourth-degree term vanishes exactly when $\alpha =-1$. For that value the fifth-degree coefficient is $-1/90\ne 0$, soAt $\alpha =1$ these are formal residual conditions only: the characteristic polynomial has a double root at one, so the method is not convergent. Canceling its common polynomial factor would change the recurrence and give a different order calculation. This distinction is relevant to part (b).

Consistency was established in part (a). By the Dahlquist equivalence theorem, convergence additionally requires the root condition for a multistep method. Since $\rho (w)=(w-1)(w-\alpha )$, its roots lie in the closed unit disk with each unit-modulus root simple exactly forThe endpoint $\alpha =-1$ is zero-stability although it has a simple parasitic root at minus one; $\alpha =1$ has a repeated root at one and is excluded.

For every parameter in this convergent range, the leading coefficient of $\sigma$ is $(5+\alpha )/12>0$, whileSince $\sigma (w)\to +\infty$ as $w\to -\infty$, it has a real root $w_{\ast }<-1$. Its other real root lies between minus one and one because $\sigma (1)=1-\alpha >0$, so the exterior root is simple.

On the Dahlquist test equation, the amplification polynomial of a multistep method is $P_z(w)=\rho (w)-z\sigma (w)$, with $z=h\lambda$. For $|z|\to \infty$, the equivalent polynomial $-P_z/z=\sigma -\rho /z$ tends coefficientwise to $\sigma$, and its leading coefficient stays nonzero. Choose a small closed disk about $w_{\ast }$ lying wholly outside the unit disk and containing no other root of $\sigma$. On its boundary $|\sigma |$ has a positive minimum, while $|\rho /z|$ is smaller for every sufficiently large $|z|$. Rouché's theorem then puts one root of $P_z$ inside this exterior disk for all such complex $z$.

That root violates the stability root condition for a multistep method. Hence there is a finite radius $R_{\alpha }$ beyond which no point belongs to the absolute stability domain:This proves boundedness in every complex direction, rather than only excluding A-stability. It is an instance of exterior roots of the derivative polynomial bound multistep stability.

For fixed real $\kappa$, consider initial data $u_0\in L^2(0,1)$ and a solution continuous into $L^2$ at time zero, with homogeneous Dirichlet boundary conditions for positive times. Existence, uniqueness and continuous dependence in this norm give Hadamard well-posedness.

To construct a solution, put $v(x,t)=e^{\kappa x}u(x,t)$. Substitution givesLet $b_j=2{\int }_0^1{e}^{\kappa x}{u}_0(x)\sin (j\pi x)dx$. The Fourier sine series solution isMultiplication by $e^{\pm \kappa x}$ is bounded on $L^2(0,1)$. Completeness of the sine basis gives convergence to $u_0$ in $L^2$ at time zero, while the exponential decay gives smoothness for positive times and the required equation and boundary values. For a classical solution continuous up to the endpoints at time zero, the usual initial-boundary compatibility conditions are additionally required.

For the energy method, integration by parts givesBoth boundary terms vanish. The same identity applies to the difference of two solutions, proving uniqueness and continuous dependence. The Poincaré inequality further giveswith the identical estimate for differences of initial data. This Dirichlet convection-diffusion contraction holds for every real $\kappa$, and its energy constant is independent of the drift strength.

Write $h=\Delta x$, impose $u_0=u_{M+1}=0$, and use the discrete inner product $⟨u,v{⟩}_h=h{\sum }_{m=1}^M{\overline{u}}_m{v}_m$. Discrete integration by parts for the second difference givesFor the centered drift, reindexing with the zero endpoints givessince the two sums are complex conjugates. This is the fact that centered discrete advection is skew-adjoint. The semidiscrete energy estimate is thereforeThe finite-dimensional linear system has a unique solution for all times, and $\Vert u(t){\Vert }_h\le \Vert u(0){\Vert }_h$, uniformly in the mesh and in real $\kappa$. Differences of two semidiscrete solutions obey the same bound, which is the required stability.

One may also use the smallest eigenvalue ${\lambda }_h=4h^{-2}{\sin }^2(\pi /[2(M+1)])$ of the negative second-difference operator to obtain $\Vert u(t){\Vert }_h\le e^{-{\lambda }_{h}t}\Vert u(0){\Vert }_h$. This conclusion concerns the method of lines; a later choice of time integrator has its own stability conditions.

For a Runge-Kutta method with real stage matrix $A=(a_{ij})$ and weights $b_i$, defineThe method has algebraic stability of a Runge-Kutta method ifThe matrix condition means $z^{T}Mz\ge 0$ for every real vector $z$. Its significance is nonlinear contractivity: an algebraically stable method, on well-defined stage solutions, is B-stable for a dissipative vector field satisfying $⟨f(t,y)-f(t,z),y-z⟩\le 0$. This is distinct from scalar A-stability, which concerns the linear test equation. The contractivity proof in part (b) explains the matrix criterion directly.

The stated hypotheses are precisely the algebraic criterion in part (a). To prove its stability content, compare two Runge-Kutta method steps for the same dissipative vector field. Let $d=y_n-z_n$ be the input difference, $D_i=Y_i-Z_i$ the stage differences, and $F_i=f(t_n+c_{i}h,Y_i)-f(t_n+c_{i}h,Z_i)$. ThenExpanding the squared norm of $d_{+}$ givesIn each term of the linear sum, replace $d$ by $D_i-h{\sum }_j{a}_{ij}{F}_j$. Symmetrizing the resulting double sum proves the Runge-Kutta contractivity identity:Dissipativity makes $⟨D_i,F_i⟩\le 0$, and the weights are nonnegative. For the second term, positive semidefiniteness permits a factorization $M=C^{T}C$, givingFor $h\ge 0$ both contributions in the distance-change identity are nonpositive. Therefore $\Vert y_{n+1}-z_{n+1}\Vert \le \Vert y_n-z_n\Vert$. This proves the contractivity assertion behind algebraic stability, usually called the Butcher contractivity theorem. It assumes the compared implicit stage equations have solutions; the matrix condition by itself should not be substituted for a separate stage-solvability argument.

For the fourth-order two-stage Gauss collocation method, $b_1=b_2=1/2$, $a_{11}=a_{22}=1/4$, $a_{12}=1/4-\sqrt{3}/6$ and $a_{21}=1/4+\sqrt{3}/6$. The weights are positive. For the matrix defining algebraic stability,andHence $M=0$ is positive semidefinite, so the method is algebraically stable. Its two-stage contractivity estimate follows from part (b). The PDF has $a_{22}=1/4$; the extra $\sqrt{3}/6$ attached to that entry in the converted TeX is a transcription error.

Put $h=\Delta x$, $k=\Delta t$ and $\mu =k/h$, consistent with $u_t=u_x$ and the exact transport $u(x,t)=F(x+t)$. The PDF's middle implicit value is $u_m^{n+1}$; the converted TeX incorrectly places it at $m+1$.

Write the new-level coefficients at offsets $-1,0,1$ as $a=\mu (1+\mu )/2$, $b=(1+\mu )(2-\mu )$ and $c=(1-\mu )(2-\mu )/2$, and the old-level coefficients as $d=2-\mu$, $e=1+\mu$. Their sums are $a+b+c=d+e=3$. To compute the local truncation error, expand the exact values $F(x_m+t_n+h(\mu +j))$ at the new level and $F(x_m+t_n+jh)$ at the old level. Equivalently the residual symbol isThe coefficients through $z^3$ vanish, whileFor example, the first three moment identities are ${\sum }_j{a}_j(\mu +j{)}^r=e$ for $r=1,2,3$, where $a_{-1}=a,a_0=b,a_1=c$; the fourth-moment difference is $\mu (\mu -2)(\mu -1)(\mu +1)$. Thus the exact residual isNormalizing the equation by $3k$ givesFor fixed nonzero $\mu \notin \{-1,1,2\}$, the method has order three. Its one-step residual is fourth degree, but the differential-equation truncation error is third order because of division by the time step.

There are exact transport exceptions. At $\mu =-1$ the equation reduces to $u_m^{n+1}=u_{m-1}^n$. At $\mu =2$ it gives $u_m^{n+1}=u_{m+2}^n$. At $\mu =1$, the exact shift $u_m^{n+1}=u_{m+1}^n$ satisfies the implicit equation, whose symbol is invertible as shown below. These are exact for every grid profile, not merely fourth-order exceptions. At $\mu =0$ the equation is the zero-time-step identity. These properties are summarized by the rational implicit advection stencil with exact shift exceptions.

For Fourier stability analysis, take spatial modes $e^{im\theta }$ and putA bounded implicit step requires $D$ to stay away from zero, and von Neumann stability analysis requires $|G|\le 1$ for every frequency. Expanding the squared moduli and using $\cos 2\theta =2{\cos }^2\theta -1$ givesConsequently, apart from denominator invertibility, the necessary and sufficient sign condition isIts real solution set is $(-\infty ,-1]\cup [0,1]\cup [2,\infty )$.

To check the denominator rather than assume it, noteIf $\mu \ne 1/2$, a zero can occur only at zero or the Nyquist frequency. The zero-frequency value is nonzero, and the Nyquist zeros occur only at $\mu =(1\pm \sqrt{3})/2$, both outside the proposed stable intervals. At $\mu =1/2$, $D=9/4+(3/4)\cos \theta \ge 3/2$, also nonzero. Continuity on the frequency circle now gives a positive lower bound on $|D|$ for each allowed $\mu$.

The implicit convolution operator is therefore boundedly invertible on the square-summable grid sequences. By the Plancherel theorem, $|G|\le 1$ givesOutside these intervals the modulus difference is negative at every nonzero frequency where $D$ is nonzero. There is an interval on which $|G|>1$, giving arbitrarily large growth of powers of the step operator, so the method is unstable; a zero of $D$ additionally makes the step ill posed.

For the physical choice $k,h>0$, the stable range is $0<\mu \le 1$ or $\mu \ge 2$; zero is the limiting zero-step case. The exact-shift values have unit gain and are included.

The necessary formulation uses a real Hilbert space $V$ of functions satisfying the essential boundary conditions, a continuous symmetric bilinear form $a$, and a continuous linear functional $\ell$. Assume constants $C,c>0$ withThe latter is coercivity, the positive-definiteness condition needed in the chosen complete space. On the domain of the strong differential operator, integration by parts and the boundary conditions identify $a(v,w)=⟨\mathcal{L}v,w⟩$; this form extends to $V$. Let $\ell (w)=⟨f,w⟩$ for admissible forcing. DefineWhen $u$ has the regularity to belong to the strong operator domain, this is $\mathcal{L}u=f$ with its boundary conditions. For complex spaces use a Hermitian form and $-2\mathrm{Re}\ell (v)$ to make the functional real.

Here is an existence proof by minimization. Continuity of $\ell$ and coercivity give $I(v)\ge c\Vert v{\Vert }_V^2-2\Vert \ell {\Vert }_{V^{\ast }}\Vert v{\Vert }_V$, so $m={\inf }_{V}I$ is finite and a minimizing sequence $v_n$ is bounded. The quadratic parallelogram identity isSince its left side is at least $m$,Coercivity makes $v_n$ Cauchy in $V$. Completeness gives $v_n\to u\in V$, and continuity of the form and forcing gives $I(u)=m$.

For every $w\in V$, differentiation along $u+tw$ yieldsThus the Euler-Lagrange equation is exactly the stated weak solution equation. Conversely, if $u$ satisfies it, expansion givesThis proves that the weak solution is the unique global minimum and that two weak solutions coincide. It proves the symmetric coercive case underlying the Lax-Milgram theorem directly, rather than taking existence as a theorem name. The framework is symmetry and coercivity in quadratic energy minimization.

The qualifications matter if positive definiteness is read more weakly. For a nonsymmetric operator the derivative of $⟨\mathcal{L}v,v⟩$ uses $(\mathcal{L}+{\mathcal{L}}^{\ast })/2$, so it need not give $\mathcal{L}u=f$. For example, on smooth periodic real functions $\mathcal{L}=1+{\partial }_x$ has $⟨\mathcal{L}v,v⟩=\Vert v{\Vert }_2^2$, but the functional's equation is $u=f$, not $u+u_x=f$. Even a symmetric operator can be strictly positive without being coercive: on l2 sequence space, $(Lv{)}_n=v_n/n$ with $f_n=1/n$ would require a solution $u_n=1$, which is not square summable. The associated functional has values $-{\sum }_{n=1}^{N}1/n$ on vectors with the first $N$ entries equal to one, so it has no minimum in that space. This is strict positivity without coercivity can fail variational solvability. Under the usual symmetric coercive meaning of positive definiteness, the proof above supplies every assertion requested.

Let $\Omega =(0,1{)}^2$ and $V=H_0^1(\Omega )$, the closure of smooth compactly supported functions in the Sobolev space $H^1$. Its elements have zero boundary trace. SetThe form is symmetric and continuous by Cauchy-Schwarz. For a smooth function with zero boundary values, $v(x,y)={\int }_0^x{v}_x(s,y)ds$ gives $\Vert v{\Vert }_{L^2}\le \Vert v_x{\Vert }_{L^2}\le \Vert \nabla v{\Vert }_{L^2}$ after integration. Density extends this Poincaré inequality to $H_0^1$. Therefore the gradient norm makes $V$ complete and equivalent to its usual $H^1$ norm, withAll the hypotheses of part (a) hold. Green's identity identifies the weak solution of the Dirichlet problem asFor smooth $u$ this is $-{\nabla }^{2}u=f$ with zero boundary values; in general it defines the unique weak solution without requiring pointwise second derivatives.

For the Ritz method, choose a finite-dimensional conforming space $V_h\subset H_0^1(\Omega )$ with basis ${\varphi }_1,\dots ,{\varphi }_N$ and put $u_h={\sum }_j{U}_j{\varphi }_j$. Minimize $I(u_h)$ in the coefficients. Its derivative in each basis direction gives the finite element method equationsThe stiffness matrix is symmetric positive definite: for a nonzero coefficient vector $U$, the corresponding basis combination is nonzero and $U^{T}KU=\Vert \nabla u_h{\Vert }_2^2>0$. Hence the Ritz system has a unique solution.

For example, triangulate the square and use continuous piecewise-affine nodal basis functions at the interior vertices, with boundary vertex values fixed to zero. On a triangle $T$ of area $A_T$, let ${\lambda }_i$ be its barycentric basis functions. If its counterclockwise vertices are $(x_i,y_i)$, define ${\beta }_i=y_j-y_k$ and ${\chi }_i=x_k-x_j$ cyclically. Then $\nabla {\lambda }_i=({\beta }_i,{\chi }_i)/(2A_T)$, giving the directly assemblable element matrices and loadsIf $f$ is constant on $T$, $F_i^T=f_T{A}_T/3$; otherwise integrate or quadrature the load. Assemble by adding element contributions with the same global vertex labels and eliminate prescribed boundary coefficients. Only overlapping basis supports interact, making $K$ sparse. This gives the Ritz-Galerkin equivalence for a symmetric coercive form in a practical finite element form.

The Euclidean matrix norm is the induced operator norm $\Vert C{\Vert }_2={\sup }_{\Vert x{\Vert }_2=1}\Vert Cx{\Vert }_2$. By the spectral theorem, a real symmetric matrix has $A=Q\Lambda Q^T$ with $Q$ orthogonal and $\Lambda =\mathrm{diag}({\lambda }_1,\dots ,{\lambda }_m)$. HenceOrthogonal transformations preserve the Euclidean norm, and the induced norm of a real diagonal matrix is the maximum absolute diagonal entry. Therefore, for $t\ge 0$,This gives the requested inequality with equality. If an exponent $\beta <\mu (A)$ were to work with prefactor one for every $t$, choose a unit eigenvector for $\mu (A)$. Its image has length $e^{t\mu (A)}>e^{t\beta }$ for every $t>0$, a contradiction. Thus $\mu (A)$ is the smallest possible exponent.

For symmetric matrices the spectral abscissa equals the Euclidean logarithmic norm. Symmetry is essential to this particular exact norm formula; a general non-normal matrix can have transient growth not controlled with prefactor one by its spectral abscissa.

For the ordered matrix exponential product $\Phi (t)=e^{tA}{e}^{tB}$, differentiation givesThe matrix commutator here is $[X,Y]=XY-YX$. Put $E(t)=\Phi (t)-e^{t(A+B)}$. Then $E(0)=0$ andMultiplying on the left by $e^{-t(A+B)}$ gives an exact product derivative. Integration from zero to $t$, followed by multiplication by $e^{t(A+B)}$, proves the exponential product defect identity:No commutation of $A$ and $B$ was assumed. The ordered product is the usual Lie-Trotter splitting; the PDF calls this product the Beam-Warming splitting. If $[A,B]=0$, the commutator term vanishes and the product is exact. In general its small-step expansion has local error $t^2[A,B]/2+O(t^3)$.

Let $a=\mu (A)$, $b=\mu (B)$ and $c=\mu (A+B)$. The sum is also symmetric negative definite. Part (a) gives the exponential norm bounds, while the triangle inequality and submultiplicativity giveApplying these estimates to the exponential product defect identity yieldsIf $c<a+b$, evaluate the scalar integral to obtain the decaying bound for symmetric negative matrix splitting:The sign condition is natural: for every unit vector $v$, $v^T(A+B)v\le a+b$, so the Rayleigh quotient gives $c\le a+b$.

When $c=a+b$, the integrand equals $e^{tc}$ independently of $x$, givingThis also follows by taking the continuous limit of the quotient as $a+b-c\to 0$. There is no division-by-zero exception in the original integral argument. All exponents are negative, so these bounds decay for large $t$; they are coarse bounds and need not reproduce the sharper quadratic small-$t$ error.