Finite Difference Methods  Example Sheet 2  A. J. Mestel 1996

1 Show that the 3-dimensional ADI method given below is an unconditionally stable approximation to the equation u_{t} = ∇^{2}u.

\begin{matrix}(1 − \frac{1}{2}rδ_{x}^{2})U_{lmn}^{∗j+1} & = [1 + r(\frac{1}{2}δ_{x}^{2} + δ_{y}^{2} + δ_{z}^{2})] U _{lmn}^{j} \\ (1 − \frac{1}{2}rδ_{y}^{2})U_{∗lmn}^{∗j+1} & = [1 + r(\frac{1}{2}δ_{x}^{2} + \frac{1}{2}δ_{y}^{2} + δ_{z}^{2})] U _{lmn}^{j} + \frac{1}{2}rδ_{x}^{2}U_{lmn}^{∗j+1} \\ (1 − \frac{1}{2}rδ_{z}^{2})U _{lmn}^{j+1} & = [1 + \frac{1}{2}r(δ_{x}^{2} + δ_{y}^{2} + δ_{z}^{2})] U _{lmn}^{j} + \frac{1}{2}rδ_{x}^{2}U_{lmn}^{∗j+1} + \frac{1}{2}rδ_{y}^{2}U_{∗lmn}^{∗j+1} \end{matrix} \} (1)

2 Show that the implicit θ-method for the hyperbolic equation u_{tt} = u_{xx} can be written

δ_{t}^{2}U_{n}^{j} = s^{2}δ_{x}^{2}(1 + θδ_{t}^{2})U_{n}^{j} ,  (2)

where s = k/h. Show that the truncation error of this scheme is O(h^{4} + k^{4}) for the particular value θ = (s^{2} − 1)/12s^{2}. Would this be a good value to choose in practice?

3 The ‘Box scheme’ for solving the advection equation u_{t} + V u_{x} = 0 (V constant) is similar to the implicit Crank-Nicolson scheme, only it treats t and x derivatives similarly:

\frac{1}{2} [\frac{(U_{n+1}^{j+1} − U_{n+1}^{j})}{k} + \frac{(U_{n}^{j+1} − U_{n}^{j})}{k}] + \frac{V}{2} [\frac{(U_{n+1}^{j+1} − U_{n}^{j+1})}{h} + \frac{(U_{n+1}^{j} − U_{n}^{j})}{h}] = 0 . (3)

Calculate the amplification factor λ for an x-dependence like e^{inξ} and deduce that the scheme is conservative. Obtain the phase shift in the form

Arg(λ) = −2 tan^{−1} [q tan(ξ/2)] , where q = \frac{kV}{h} .  (4)

Observe that the low harmonics (ξ small) are modelled well by the finite difference approximation. Show also that when the Courant number q = 1, the phase factor is exactly correct for all ξ. Deduce that the approximation is exact in this case.

4 Examine the stability of the two-step Lax-Wendroff method for the conservation equation \underline{u}_{t} + [\underline{f}(\underline{u})]_{x} = 0, namely

\begin{matrix}\underline{U} _{n+1/2}^{j+1/2} & = \frac{1}{2} [\underline{U} _{n}^{j} + \underline{U} _{n+1}^{j}] − \frac{1}{2}s [\underline{f} _{n+1}^{j} − \underline{f} _{n}^{j}] \\ \underline{U} _{n}^{j+1} & = \underline{U} _{n}^{j} − s [\underline{f} _{n+1/2}^{j+1/2} − \underline{f} _{n−1/2}^{j+1/2}] \end{matrix} \}  (5)

for the case when \underline{f}(\underline{u}) = A\underline{u}, where A is a constant matrix. Writing \underline{U}_{n}^{j} = (λ)^{j}e^{inξ}\underline{V}, show that

[2s^{2} sin^{2}(\frac{1}{2}ξ)A^{2} + is sin(ξ)A + (λ − 1)I] \underline{V} = 0 .  (6)

Assuming the eigenvectors of A are complete, the only non-zero solutions to (6) occur when \underline{V} is an eigenvector of A, so that

2s^{2} sin^{2}(\frac{1}{2}ξ)μ^{2} + is sin(ξ)μ + (λ − 1) = 0 ,  (7)

where μ is an eigenvalue of A. Assuming all values of μ are real, show that the stability condition is s max | μ |≤ 1. If, conversely, μ takes complex values (in conjugate pairs, of course), show that the low harmonics (ξ → 0) grow exponentially for any value of s. Interpret this behaviour in terms of the existence of characteristics and the underlying nature of the partial differential equation.

5 Show that for a square grid with steplength h, the Laplacian may be approximated by

\frac{1}{2h^{2}} [u_{m+1,n+1} + u_{m+1,n−1} + u_{m−1,n−1} + u_{m−1,n+1} − 4u_{m,n}] = (∇^{2}u)_{m,n} + O(h^{2}) . (8)

Examine the truncation errors of (8) and the usual 5-point formula (∇^{2} ≈ (δ_{x}^{2} + δ_{y}^{2})/h^{2}) and hence find a linear combination of the two which enables an O(h^{4}) approximation to ∇^{2}u = f to be made.

6 Consider the solution of the Poisson equation with Dirichlet boundary conditions,

u_{xx} + u_{yy} = f  in D,  u = g  on a curved boundary C.  (9)

Using a square mesh, we need to approximate the Laplacian near the boundary. Suppose that the curved boundary intersects the grid at the points A and B whose coordinates are (x, y) = (M + θ, N)h, and (x, y) = (M, N + ϕ)h respectively, where 0 < θ ≤ 1, 0 < ϕ ≤ 1, as in the figure. Calculate approximations to u_{xx} and u_{yy} at the mesh point (Mh, Nh) and hence show that a reasonable approximation to (9) at that point is

\frac{2g(A)}{θ(1 + θ)} + \frac{2g(B)}{ϕ(1 + ϕ)} + \frac{2U_{M−1,N}}{1 + θ} + \frac{2U_{M,N−1}}{1 + ϕ} − 2 (\frac{1}{θ} + \frac{1}{ϕ}) U_{M,N} = h^{2}f_{M,N} , (10)

where g(A) denotes g evaluated at the point A. What is the order of the truncation error of (10)?