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

1. Given a function u(x) and writing u_{n} = u(nh), show that

\frac{12u_{1/3} − 3u_{2/3} − 9u_{0}}{2h} = u′(0) + O(h^{2}) .

2. The steady equations of motion of inviscid, isentropic compressible fluid flow are

\begin{matrix}ρ(u · ∇)u + ∇p & = 0 \\ u · ∇ρ + ρ∇ · u & = 0 \end{matrix} \}  (1)

where u is the velocity vector and the pressure, p, is a function of the density, ρ, only: p = p(ρ). Defining the local sound speed c by c^{2} = \frac{dp}{dρ}, show that these equations imply that [skip this bit if you have no knowledge of fluid dynamics]

u · ∇(\frac{1}{2}u · u) = c^{2}∇ · u .  (2)

If the flow is both two-dimensional and irrotational (∇ ∧ u = 0), then u can be written in terms of a potential, ϕ(x, y), so that u = (ϕ_{x}, ϕ_{y}, 0). Obtain an equation for ϕ, and determine the conditions for it to be hyperbolic, elliptic or parabolic, commenting on the physical meaning of the results.

3. The equation u_{t} = u_{xx} may be modelled by a two-step method which is centred in both space and time, namely

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

What is the order of the truncation error of this scheme?

Use the Fourier method to investigate the stability of (3) for all positive values of r = k/h^{2}.

4. The stability behaviour of the method (3) can be improved by replacing the term 2U_{n}^{j} by the average (U_{n}^{j+1} + U_{n}^{j−1}):

U_{n}^{j+1} − U_{n}^{j−1} = 2r (U_{n+1}^{j} − U_{n}^{j+1} − U_{n}^{j−1} + U_{n−1}^{j}) .  (4)

Calculate the truncation error of this scheme up to and including terms of O(k^{2}) and O(h^{2}).

Under what circumstances is this approximation consistent? i.e. as h → 0 and k → 0, when does the truncation error of the finite difference approximation tend to zero?

5. Hard: Use the Matrix Method to examine the stability of the scheme

(1 + θ) [\frac{U_{n}^{j+1} − U_{n}^{j}}{k} ] − θ [\frac{U_{n}^{j} − U_{n}^{j−1}}{k} ] = \frac{U_{n+1}^{j} − 2U_{n}^{j} + U_{n−1}^{j}}{h^{2}}  (5)

where 0 ≤ θ ≤ 1. (Note θ = −\frac{1}{2} corresponds to (3).)

6. Pollutant is introduced at a rate Q(x, y) into a river flowing with velocity u(x) in the y-direction. Equilibrium between downstream advection and cross-channel diffusion is reached when the pollutant concentration C(x, y) satisfies

uC_{y} = DC_{xx} + Q ,  where D is constant.  (6)

The river occupies −d < x < d, and a parabolic velocity profile u = u_{0}(d^{2}−x^{2}) is assumed. It is desired to solve this equation subject to boundary conditions on x = ±d and an initial condition at y = 0 by an explicit method. If the steplength used across the river is 1/20 of its width, what is the maximum steplength along the river which will guarantee stability of the algorithm?

Should the direction of solution be upstream or downstream? Is the method a good one?

7. The three-dimensional diffusion equation

u_{t} = ∇^{2}u = u_{xx} + u_{yy} + u_{zz}  (7)

is to be solved on a rectangular grid with step-lengths (h_{x}, h_{y}, h_{z}) using the θ-method

U_{lmn}^{j+1} − U_{lmn}^{j} = k [\frac{1}{h_{x}^{2}} δ_{x}^{2} + \frac{1}{h_{y}^{2}} δ_{y}^{2} + \frac{1}{h_{z}^{2}} δ_{z}^{2}] (θU_{lmn}^{j+1} + (1 − θ)U_{lmn}^{j}) ,  (8)

in the usual notation. Use the Fourier method to show that this scheme is unconditionally stable for 1 ⩾ θ ⩾ \frac{1}{2} while if 0 ⩽ θ < \frac{1}{2} it is stable provided

k \{\frac{1}{h_{x}^{2}} + \frac{1}{h_{y}^{2}} + \frac{1}{h_{z}^{2}} \} ⩽ \frac{1}{2(1 − 2θ)} .  (9)