Past exam of the numerical solution of differential equations course of the University of Oxford / 2011 / A

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