Past exam of the numerical solution of differential equations course of the University of Oxford / 2010 (B21, DMAB 2639, DCOM 2639, DMCN 2639, DMST 2639)

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$.

The coefficient signs alone do not imply the assertion as printed. Two interior nodes can have rows $L_h{\psi }_1={\psi }_2-{\psi }_1$ and $L_h{\psi }_2={\psi }_1-{\psi }_2$, with a separate boundary node. Taking interior values one and boundary value zero satisfies every stated coefficient inequality and $L_h\psi =0$, but its maximum is not on the boundary. Connectedness of the geometric domain does not ensure connectedness of this stencil graph.

The usual extra assumption is that every interior node has a path of positive stencil coefficients to a prescribed boundary node. Under this assumption let $M\ge 0$ be the global maximum and suppose an interior node $p$ attains it. ThenAll summands are nonpositive, so equality forces ${\psi }_r=M$ at every neighbor with $c_{p,r}>0$. Propagate along a path to the boundary to find a boundary node with value $M$. This proves the discrete maximum principle under boundary-reaching stencil connectivity, and exhibits why the literal hypotheses need that qualification.

A diameter-only bound also needs the normalization of a Laplacian stencil, beyond its coefficient signs. For example, a single interior row $L_h{e}_p=\epsilon (e_a-e_p)$ with $e_a=0$ allows $|L_h{e}_p|=1$ and $|e_p|=1/\epsilon$ on a domain of fixed diameter. Thus arbitrary rescaling of the printed coefficients rules out a universal diameter-only estimate.

For the intended normalized Laplacian stencil, assume $L_{h}1=0$ and $L_h\varphi =4$ for ${\varphi }_p=|x_p-x_0{|}^2$, with $x_0\in \Omega$, as holds for the standard two-dimensional central stencil. Let $D=\mathrm{diam}\Omega$. Then $0\le {\varphi }_p\le D^2$. Use the truncation error convention $T_p=L_h{u}_p-f_p$, so $L_h{e}_p=T_p$; reversing the convention does not change the absolute-value argument.

The two functions ${\psi }_p^{\pm }=\pm e_p+(T/4){\varphi }_p$ satisfy $L_h{\psi }^{\pm }\ge -T+T=0$. On the prescribed boundary $e_a=0$, so their maximum boundary values are at most $T{D}^2/4$. Apply the discrete maximum principle from part (a), with its boundary-connectivity assumption. Since $\varphi \ge 0$,This is the standard quadratic-barrier error estimate for a discrete Laplacian. The calculation identifies exactly the additional normalization needed for the bound intended by the question.

The displayed estimate requires a nonnegative barrier. As printed, a counterexample is obtained with two interior nodes, each having row $L_h{v}_p=v_a-v_p$ to one boundary node. Set $e_a=0$, $e_p=1$, and ${\varphi }_a=-1$, ${\varphi }_p=-2$. Both interior rows have $|T_p|=1$ and $L_h{\varphi }_p=1$. Put one node in each of $J_1,J_2$ and take $T=S=d_1=d_2=1$. Then the proposed right-hand side is $-1$ while $E=1$.

For the intended use, suppose ${\varphi }_p\ge 0$ and the stencil has the maximum principle established in part (a). Set $C=\max \{T/d_1,S/d_2\}$. On either subregion,At boundary nodes these functions equal $C{\varphi }_a$. Hencewhich proves the intended discrete elliptic barrier estimateFor a normalized stencil with $L_{h}1=0$ and a sign-changing barrier, shifting it by its minimum gives the valid general replacement $E\le C({\max }_a{\varphi }_a-{\min }_p{\varphi }_p)$. The supplied barrier in part (d) is positive, so the printed estimate applies there without a sign problem.

The given stencil has nonnegative neighbor weights, zero row sums and a path from every unknown node to the Dirichlet boundary. Thus the discrete maximum principle applies, including at the reflected Neumann edge. The barrier is positive throughout the square.

At an ordinary interior node, central differences differentiate this quadratic exactly, giving $L_h\varphi =2+2=4$. At $x=1$, the reflected horizontal part isand the vertical part is two. Thus $L_h\varphi =4+4/h$ there. Part (c) therefore permitsOn the Dirichlet boundary, the largest barrier value is $4+1/4=17/4$. ConsequentlyThe $O(h)$ truncation error at the Neumann edge is divided by the barrier response $O(h^{-1})$, producing a contribution of order $h^2$. Thus its lower local order does not reduce this global maximum-norm error bound.

With zero boundary values eliminated, the diagonal of $L_h$ is $-4I/h^2$. The Jacobi method is thereforewhich averages the four neighboring components. The sine addition identity shows that the supplied product vectors are eigenvectors, withThe discrete sine vectors form a basis, so these are all eigenvalues. Their largest modulus isattained at the lowest and highest sine frequencies. Thus the asymptotic error reduction factor per iteration is $\cos (\pi h)$. If the convergence rate is expressed logarithmically, $R=-\log \rho (J)\sim {\pi }^2{h}^2/2$ as $h\to 0$. A fixed error reduction therefore takes $O(h^{-2})$ Jacobi sweeps. For the special grid $m=2$, $J=0$ and one sweep solves its single interior equation.

The endpoint $\alpha =0$ in the printed nonnegative-parameter hypothesis must be excluded. We first establish the intended result for $\alpha >0$, then describe that endpoint exactly.

Write $X,Y$ for the two directional matrices; they commute because they act on separate tensor-product coordinates. At a limit $V$, the intermediate vector also has a limit $\tilde{V}$. Subtract the two limiting equations in the form that cancels $X\tilde{V}$ and $YV$; this gives $2\alpha (\tilde{V}-V)=0$. Hence $\tilde{V}=V$ for $\alpha >0$. Substitution then gives $(X+Y)V=h^{2}f$, precisely the required discrete equation.

Eliminating the intermediate vector gives the alternating direction implicit iteration matrixThe sine-product basis simultaneously diagonalizes $X,Y$, with eigenvalues $-{\xi }_p$ and $-{\eta }_q$, whereThusFor $\alpha >0$, each factor has modulus strictly less than one. Hence $\rho (G)<1$, proving convergence from every starting vector to the unique discrete solution.

At $\alpha =0$ the two factors are both $-1$, so $G=I$. The forcing terms cancel in the composed iteration as well, and $V^{k+1}=V^k$ for every starting vector. An incorrect initial vector therefore remains incorrect forever, even though the sequence is constant and has a limit. This disproves both the literal fixed-point implication and the claimed general convergence-to-solution assertion at zero; strict positivity of $\alpha$ is the necessary repair.