Past exam of the numerical solution of differential equations course of the University of Oxford

Original solutions to both sections and all six questions of the Numerical Solution of Differential Equations paper sat on 22 June 2011. Source: released Mathematics papers, Part B 2007 11.pdf. The real-triangularization assertion in question 6 is corrected to complex triangularization. Statements are not reproduced.

Use the normalized local truncation errorThe Taylor theorem expansions about $x_n$ areSubstitution gives exactlyThus the theta method has order two exactly for $\theta =1/2$, and otherwise order one, as a statement about the method on a sufficiently smooth class of problems. The leading term vanishes at one half, and a function with nonzero third derivative shows that the next term does not vanish universally. As usual, a uniform order estimate assumes bounded relevant derivatives; mere existence of $y^{\prime \prime \prime }$ on the open interval gives the stated pointwise remainders but not, by itself, a uniform bound. Special exact solutions can of course have zero truncation error even for a first-order method.

Put $e_n=y(x_n)-y_n$. Subtracting the numerical update from the exact update with its local truncation error givesThe Lipschitz condition in the solution variable and the triangle inequality implyThe assumed positive denominator permits division. With $K=[1+L(1-\theta )h]/(1-L\theta h)$ and $e_0=0$, repeated substitution of this Lipschitz one-step error bound yieldsHere $K>1$ because $L,h>0$, so the displayed geometric series quotient is well defined. For implicit steps, the same condition makes the scalar update map a contraction mapping, ensuring a unique numerical next value.

On the Dahlquist test equation, set $z=h\lambda$. The theta method has stability functionWith the decaying-solution convention, its region of absolute stability is $|R(z)|<1$. Writing $z=x+iy$, squaring numerator and denominator givesFor $0\le \theta <1/2$ this is the open disk with centre $-1/(1-2\theta )$ and radius $1/(1-2\theta )$. At $\theta =1/2$ it is the open left half-plane. For $\theta >1/2$ it is the exterior of the closed disk with centre $1/(2\theta -1)$ and radius $1/(2\theta -1)$. The pole $z=1/\theta$, when present, lies outside the stability region.

Every $z$ with $\mathrm{Re}z<0$ satisfies the inequality exactly for $\theta \ge 1/2$. If $\theta <1/2$, sufficiently large negative real $z$ fails it. ThereforeIf absolute stability is instead defined by bounded amplification, the boundary $|R|=1$ is included; this convention does not change the A-stability conclusion for the open left half-plane.

In the normalization used by the supplied coefficient expansion, the local truncation error of the linear multistep method isThe nonzero sum in the denominator is an explicit hypothesis. Consistency of a numerical method requires $T_n\to 0$ uniformly on a fixed finite time interval for every sufficiently smooth exact solution as $h\to 0$. Equivalently, the two polynomial identities areThese follow by inserting constant and affine functions, or by setting the first two Taylor coefficients $C_0,C_1$ to zero.

Zero-stability means that perturbations of the starting values of the zero-step recurrence remain bounded by a fixed multiple of their initial size, uniformly as the number of steps tends to infinity on a fixed time interval. In particular solutions of ${\sum }_{j=0}^k{\alpha }_j{e}_{n+j}=0$ obey ${\max }_{n\le N}|e_n|\le C{\max }_{j<k}|e_j|$ with $C$ independent of $N$.

The equivalent root condition for a multistep method isA root outside the unit circle gives exponential growth of a recurrence relation mode; a repeated root on it gives polynomial growth. Repeated roots strictly inside it are harmless because their polynomial factor is dominated by decay.

The first characteristic polynomial factors asThe two quadratic roots have product one. Both can lie in the closed unit disk only if both have modulus one, requiring $|\alpha +1|\le 2$. At $\alpha =-3$ the root one is triple and $\beta =0$; at $\alpha =1$ the root minus one is double. Both fail the root condition for a multistep method. In the intervening open interval the two quadratic roots are distinct conjugate unit-modulus roots, different from one. ThusUsing $2\beta =\alpha +3$, the supplied expansion givesThroughout the zero-stable range, $C_3\ne 0$, so every such method has order two. There is a unique higher-order choice:Here $C_3=C_4=0$ and $C_5=1/10\ne 0$, giving exactly order four, with normalized leading truncation error $h^4{y}^{(5)}/120$. But the extra roots are $-5\pm 2\sqrt{6}$, one of modulus greater than one. By the Dahlquist equivalence theorem this consistent but non-zero-stable method is not convergent for arbitrary consistent starting values.

Let ${\delta }_{xx}{U}_{ij}=(U_{i+1,j}-2U_{ij}+U_{i-1,j})/(\Delta x{)}^2$, and define ${\delta }_{yy}$ analogously. The Crank-Nicolson diffusion scheme isUse $U_{ij}^0=u_0(x_i,y_j)$ and zero values on every boundary node for every time level. This is the theta method with $\theta =1/2$ applied to the spatial finite difference method for the heat equation.

Write $s={\mu }_x+{\mu }_y$. The update can be rearranged asWhen $s\le 1$, the right side is a convex combination of old values: its coefficients are nonnegative and sum to one. Assume all old values lie in $[U_{\min },U_{\max }]$. A new maximum $v>U_{\max }$ is interior because boundary values are zero and $U_{\max }\ge 0$. At its node the left side is at least $(1+s)v-sv=v$, while the right side is at most $U_{\max }$, a contradiction. A new minimum $v<U_{\min }$ similarly makes the left side at most $v$ and the right side at least $U_{\min }$, again a contradiction.

The initial values satisfy the bounds by construction; mathematical induction provesThe new-time matrix is a strictly diagonally dominant matrix with positive diagonal and nonpositive off-diagonal entries, so the linear update is uniquely solvable. This Crank-Nicolson maximum principle under a step restriction concerns preservation of value bounds, which is a stronger requirement than scalar A-stability.

For nonnegative initial data, $U_{\min }=0$. The proved Crank-Nicolson maximum principle under a step restriction immediately gives $\boxed{U_{ij}^m\ge 0}$ for all indices. For nonpositive initial data, $U_{\max }=0$ and the same bound gives $\boxed{U_{ij}^m\le 0}$. Equivalently, apply positivity to $-U$, since the finite difference method is linear and has homogeneous Dirichlet boundary conditions. Both conclusions use ${\mu }_x+{\mu }_y\le 1$.

For the second-order boundary value problem, introduce an unknown initial slope $s$ and solve the initial value problem $u^{\prime }=v$, $v^{\prime }=g(x,u,v)$ with $u(a)=\alpha$, $v(a)=s$. Its endpoint mismatch is $F(s)=u(b;s)-\beta$. The shooting method seeks a root $F(s)=0$; in general the dependence of the flow on $s$ is nonlinear. A sign-changing bracket for a continuous $F$ allows the bisection method, provided the intervening initial value solutions exist up to $b$.

For Newton's method, solve the parameter sensitivity equation alongside the trajectory. With $z={\partial }_{s}u$ and $w={\partial }_{s}v$,Then $F^{\prime }(s)=z(b)$ and the update is $\boxed{s_{\mathrm{new}}=s-F(s)/z(b)}$ whenever $z(b)\ne 0$.

For the third-order problem, shooting from $a$ leaves both $u^{\prime }(a)$ and $u^{\prime \prime }(a)$ unknown, requiring two shooting parameters to satisfy two right-endpoint conditions. Shooting backwards from $b$ already knows $u(b)=\beta$ and $u^{\prime }(b)=\gamma$; only $u^{\prime \prime }(b)$ is unknown, with the single residual $u(a)-\alpha$. Backward shooting reduces the nonlinear solve from two scalar unknowns to one. This parameter-count advantage does not by itself guarantee numerical stability of backward integration.

The four-point midpoint third derivative uses the four values at $x_{j-1},x_j,x_{j+1},x_{j+2}$, centred at $x_{j+1/2}$. The Taylor theorem about that midpoint givesThus, for $n\ge 2$, take the $n-1$ interior equationsImpose the two value conditions exactly, $U_0=\alpha$, $U_{n+1}=\beta$, and approximate the remaining condition with the second-order backward finite difference:This gives $n$ equations for $U_1,\dots ,U_n$. The endpoint derivative formula equals $u^{\prime }(b)-h^2{u}^{\prime \prime \prime }(b)/3+O(h^3)$, so its residual tends to zero along with those of the four-point equations. Every stencil lies within $[a,b]$; no unprescribed ghost value is introduced. The construction is therefore a consistent finite difference method for all three boundary conditions.

The five-point equations for the Poisson equation arewith prescribed boundary values. If $L_h\varphi \le 0$, an interior maximal value is at most the average of its four neighbours. Since no neighbour exceeds the maximum, equality forces all four to equal it. Repeating along adjacent nodes propagates that value to the boundary. Thus the discrete maximum principle gives $\max \varphi \le {\max }_{\partial {\Omega }_h}\varphi$, in particular the usual positive-maximum bound with zero included.

For an $O(h^2)$ error estimate assume the exact solution has bounded continuous fourth partial derivatives up to the boundary. The local truncation error $\tau =L_{h}u-f$ then satisfies $\Vert \tau {\Vert }_{\infty }\le C{h}^2$ by the Taylor theorem. Put $e=u{|}_{{\Omega }_h}-U$, so $L_{h}e=\tau$ and $e=0$ on boundary nodes. The Discrete Poisson barrier on a squaresatisfies $L_{h}B=1$ exactly, $B\ge 0$ on the grid boundary, and $B\le 1/8$. Set $\epsilon =\Vert \tau {\Vert }_{\infty }$. Then $L_h(\pm e-\epsilon B)\le 0$, with nonpositive boundary values. Applying the discrete maximum principle twice givesFor nonzero smooth Dirichlet boundary data $g$, impose $U=g$ at the boundary nodes. The error still vanishes there, so the same argument and order apply, under the same fourth-derivative regularity of the exact solution. Smooth boundary values alone do not remove every possible corner-regularity issue.

Replace $L_h$ by $\boxed{{\tilde{L}}_h=L_h+I}$. At a positive interior maximum $v>0$, the neighbour average is no larger than $v$, so $({\tilde{L}}_h\varphi {)}_{jk}\ge v>0$. Thus ${\tilde{L}}_h\varphi \le 0$ cannot have such an interior maximum. The positive-maximum form of the discrete maximum principle therefore remains valid:With $e=u{|}_{{\Omega }_h}-U$ the error equation is ${\tilde{L}}_{h}e=\tau$, where $\tau =O(h^2)$ under the preceding regularity assumptions. Since ${\tilde{L}}_{h}B=1+B\ge 1$ for the same Discrete Poisson barrier on a square, the comparison functions $\pm e-\Vert \tau {\Vert }_{\infty }B$ proveThere is no loss of uniqueness: apply the positive-maximum argument to both signs of a homogeneous discrete solution. The continuous weak formulation also has the coercive bilinear form ${\int }_{\Omega }(|\nabla v{|}^2+v^2)$, so the added positive reaction term strengthens its positive energy.

Let $x_{\ast }=A^{-1}b$, $e_k=x_k-x_{\ast }$ and $G=M^{-1}N$. Subtract the exact fixed-point equation to obtain $e_k=G^k{e}_0$. Hence convergence for every initial vector is equivalent to $G^k\to 0$. If this holds, it holds after complexification; applying it to a complex eigenvector gives ${\lambda }^{k}v\to 0$, so $|\lambda |<1$. Conversely, the complex Jordan normal form expresses each power as a finite sum of terms $\left(\genfrac{}{}{0ex}{}{k}{r}\right){\lambda }^{k-r}$ from its blocks. Every term tends to zero when all eigenvalues have modulus below one. ThusThe convergence criterion for a stationary iteration needs complex triangularization or Jordan normal form. The paper's assertion that every real matrix is similar to a real upper-triangular matrix is false when nonreal eigenvalues occur; complex triangularization supplies the intended argument. In the Gauss-Seidel method, with $A=D+L+U$ split into diagonal, strictly lower and strictly upper parts, $\boxed{M=D+L,N=-U}$.

For $n\ge 4$, eliminating the ghost values gives the symmetric pentadiagonal matrixAll other off-diagonal entries vanish. The first two equations areand the last two areFor $n=3$ the matrix is $\left(\begin{array}{ccc}7& -4& 1\\ -4& 6& -4\\ 1& -4& 7\end{array}\right)$, for $n=2$ it is $\left(\begin{array}{cc}7& -4\\ -4& 7\end{array}\right)$, and for $n=1$ it is $(8)$ because both ghost contributions meet the same unknown. The diagonal entries are positive in every case, so the lower triangular matrix $D+L$ is invertible and the Gauss-Seidel iteration is well defined.

An eigenvalue $\lambda$ of $G=-(D+L{)}^{-1}U$ makes $H(\lambda )=\lambda (D+L)+U$ singular. Let ${\ell }_i$ and $r_i$ be the sums of absolute lower and upper entries in row $i$. The row's diagonal-dominance margin isInterior full rows have $d_i=6$, ${\ell }_i=r_i=5$, giving margin $|\lambda |-5$. Boundary rows have larger margins; for instance the first gives $7|\lambda |-5$, and the second gives at least $2|\lambda |-5$. Therefore $|\lambda |>5$ would make $H(\lambda )$ a strictly diagonally dominant matrix, hence nonsingular, a contradiction.

At $|\lambda |=5$ it is weakly diagonally dominant with strict boundary rows. Its nonzero nearest-neighbour entries make it an irreducibly diagonally dominant matrix. The nonsingularity theorem for that class again contradicts singularity. For the small matrices above all rows are already strict at this value. HenceThere is also a stronger conclusion. For $T=\mathrm{tridiag}(-1,2,-1)$,including the $n=1$ overlap. Since $T$ is positive definite, $A$ is a symmetric positive-definite matrix. The Gauss-Seidel convergence from diagonal dominance and positive definiteness theorem consequently gives the stronger bound $\rho (G)<1$.

Source: released examination bundle, Part B 2007 11.pdf, pages 248–255. All six questions in both halves are covered. No mark allocations are printed.

Substitute a smooth exact solution into the linear multistep method and expand about $x_n$. Its residual isConsistency of a numerical method requires the coefficients of $y$ and $h{y}^{\prime }$ to vanish. Therefore $1+a+b=0$ and $3+a=1$, giving the unique choiceFor this choice the next coefficient is $9/2-1-2=3/2\ne 0$, so the residual is $\frac{3}{2}{h}^2{y}^{\prime \prime }+O(h^3)$. The normalized local truncation error, obtained by dividing the residual by $h$, is $O(h)$. Hence the method has order one.

A zero-stable has a uniformly bounded response to small perturbations of its starting values and accumulated residuals on a fixed time interval as $h\to 0$. Equivalently, its first characteristic polynomial satisfies the root condition for a multistep method: every root has modulus at most one, and each root on the unit circle is simple.

For the consistent parameters, the recurrence for $y^{\prime }=0$ is $y_{n+3}-2y_{n+1}+y_n=0$, with characteristic polynomialLet $r=(-1-\sqrt{5})/2$. Choose $y_1=r$ and $y_2=r^2$; together with $y_0=1$ these generate exactly $y_n=r^n$. Since $|r|=(1+\sqrt{5})/2>3/2$,This parasitic mode violates the root condition and grows without a bound independent of the number of steps. Scaling its starting amplitudes gives the same unbounded amplification for arbitrarily small initial perturbations. Thus the only consistent method in the family is not zero-stable.

A method has numerical convergence if it approximates the exact solution uniformly at mesh points on every fixed finite interval as $h\to 0$, whenever its starting approximations converge to the corresponding exact starting values, for sufficiently smooth solutions of a well-posed initial-value problem.

The Dahlquist equivalence theorem states that a fixed-coefficient linear multistep method is convergent if and only if it is consistent and zero-stable, under the usual Lipschitz condition on the differential equation and convergent starting data. In this family numerical consistency forces the parameters found in part (a), and part (b) proves that choice is not zero-stable. All other parameter choices are inconsistent. Therefore no member of the family is convergent.

With $d=\Delta x$, define the semidiscrete Fourier transform byFor finite-support sequences, expand the squared modulus and use orthogonality:ConsequentlyFor a general square-summable sequence, truncate it to $|j|\le N$. The identity makes the transformed truncations Cauchy in $L^2$, and completeness defines the transform as their $L^2$ limit. Passing to the limit proves the Parseval identityThe infinite series need not converge absolutely pointwise; the $L^2$ interpretation is sufficient for the stability analysis.

Transforming the spatial shifts and putting $\xi =k\Delta x$ gives the real amplification factorFor $0\le \theta \le 1$ the denominator is at least one, so the implicit step is a bounded Fourier multiplier on ${\ell }^2$. Set $v={\sin }^2(\xi /2)\in [0,1]$ and use ${\sin }^2\xi =4v(1-v)$. The numerator is no greater than the denominator, so $G\le 1$. The other inequality isIf $\theta \ge 1/2$, the bracket is nonpositive and this always holds. If $0\le \theta <1/2$, the concave quadratic has its maximum at $v=(1-2\theta )/[2(1-\theta )]$, with value $(1-2\theta {)}^2/[4(1-\theta )]$. Therefore $|G|\le 1$ precisely under the stated conditionand without a mesh-ratio restriction when $1/2\le \theta \le 1$. The Parseval identity now givesIteration proves stability with constant one. It is conditional for $\theta <1/2$ and unconditional for $\theta \ge 1/2$. The sampled initial data have finite support, hence belong to ${\ell }^2$. This is the Fourier stability threshold for a mixed-width heat scheme.

On nonnegative states the characteristic speed is $f^{\prime }(u)=1+p{u}^{p-1}>0$, so upwind finite difference scheme uses the left spatial difference. With $\lambda =\Delta t/\Delta x$, the conservative upwind finite difference scheme isExplicitly the flux difference is $(U_j^m-U_{j-1}^m)+[(U_j^m{)}^p-(U_{j-1}^m{)}^p]$. Preservation of nonnegativity under the given Courant restriction, proved next, justifies the same left-upwind direction at every step.

Write $M_0={\sup }_x{u}_0(x)<\infty$ and assume $\lambda (1+p{M}_0^{p-1})\le 1$. Suprema are understood where the printed maxima on the infinite domain are not attained. We prove the invariant interval $0\le U_j^m\le M_0$ simultaneously with the stated speed bound.

Suppose the interval property holds at time $m$. By the mean value theorem, the divided differenceis between $1$ and $1+p{M}_0^{p-1}$; when the two states coincide define it as $f^{\prime }(U_j^m)$. The update is the convex combinationBoth coefficients are nonnegative and sum to one, so the same interval property holds at time $m+1$. It is true initially and hence for all times by induction. ThereforeThe propagated nonnegativity is essential; it keeps the flux increasing on the entire range of the numerical solution.

The convex-combination formula proved in part (i) gives $0\le U_j^{m+1}\le {\sup }_{\ell }{U}_{\ell }^m$ for every $j$. Taking the supremum and using nonnegativity givesThis is the discrete maximum principle for the nonlinear upwind step and proves the requested conditional ${\ell }^{\infty }$ stability under the propagated Courant–Friedrichs–Lewy condition.

Introduce $v=u^{\prime }$ and solve the initial-value systemThe unknown initial slope $s$ is the shooting method parameter. Integrate to one and solve the scalar residual equation $R(s)=u(1;s)-\beta =0$, using a secant iteration or a bracketed root finder when suitable data are available.

For a Newton update, integrate the sensitivity equations along the current shot:Then $R^{\prime }(s)=w(1)$ and the update is $s\leftarrow s-R(s)/w(1)$ when this derivative is nonzero. After a root has been found, the corresponding shot satisfies both endpoint conditions. This formulation does not presume that a nonlinear boundary-value problem necessarily has a unique solution.

At $x_r=rh$, the central finite difference equations arewith $U_0=\alpha$, $U_m=\beta$. For $f=f(x)$ this is a linear system $A_h{U}_h=b_h$, with row coefficientsThe right-hand side starts from $f(x_r)$ and includes $-a_1\alpha$ in its first entry and $-c_{m-1}\beta$ in its last entry. Under $h\max p<2$, both off-diagonal coefficients are negative, so $|a_r|+|c_r|=2h^{-2}<d_r$. Omitting nonexistent endpoint columns can only reduce this sum. Therefore $A_h$ is strictly diagonally dominant. In particular the linear discrete system is nonsingular.

Use the fine matrix $A_h$ and right-hand side from part (b), with diagonal $D_h$. Starting from an iterate $v$, take $\nu$ weighted Jacobi smoothing stepsForm the residual $r_h=b_h-A_{h}v$. For coarse interior indices $j=1,\dots ,m/2-1$, use full-weighting restrictionThe prolongation operator matrix $P$ is linear interpolation of a coarse correction: $(Pe{)}_{2j}=e_j$ and $(Pe{)}_{2j+1}=(e_j+e_{j+1})/2$, with zero endpoint corrections. Thus $R=P^T/2$. Choose the Galerkin coarse-grid operator $A_c=R{A}_{h}P$, solve $A_c{e}_c=R{r}_h$ exactly, and correct $v\leftarrow v+P{e}_c$. The entries of $A_h$, $R$ and $P$ above specify every matrix in the cycle. A rediscretized coarse operator is another standard choice if used consistently.

For a nonsingular coarse operator, the fine-grid iteration matrix of this pre-smoothed two-grid cycle isThe purpose of under-relaxation is to damp fine-grid oscillations, leaving smoother error for coarse-grid correction. For the diffusion-dominated one-dimensional operator, $\omega =2/3$ is the usual effective choice; the high-frequency symbol $1-\omega (1-\cos \xi )$ then has modulus at most $1/3$ for $\pi /2\le |\xi |\le \pi$. Coefficients and convection can change the best parameter, so this is a motivated choice rather than a universal optimum. Under the diagonal-dominance condition, the unweighted Jacobi matrix already has norm less than one in the maximum norm, and its convex combination with the identity remains convergent for $0<\omega <1$.