Past exam of the numerical solution of differential equations course of the University of Oxford / 2010 / B (Numerical Solution of Differential Equations II)

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.