11. The Lax-Wendroff method and conservation equations

One explicit method which attempts to get the second order terms correct is due to Lax and Wendroff. We know that

\frac{u_{n}^{j+1} − u_{n}^{j}}{k} = [u_{t} + \frac{1}{2}ku_{tt}]_{n}^{j} + O(k^{2}) .

So if u_{t} = −cu_{x} with c constant, then u_{tt} = −cu_{xt} = c^{2}u_{xx}, and we may write

\frac{u_{n}^{j+1} − u_{n}^{j}}{k} = −c\frac{Δu_{n}^{j}}{2h} + \frac{kc^{2}}{2h^{2}} δ^{2}u_{n}^{j} + O(k^{2}, h^{2})

This leads to the centred scheme (10.7F) on the last sheet

U_{n}^{j+1} = U_{n}^{j} − \frac{1}{2}qΔU_{n}^{j} + \frac{1}{2}q^{2}δ^{2}U_{n}^{j}  (11.1)

Note that this formula requires c to be constant. However, a similar result can be derived for the conservation equation u_{t} + [f(u, x)]_{x} = 0. In general, an equation of the form

\frac{∂ρ}{∂t} + ∇ · \{b\} = Q

can be interpreted as a conservation law for a physical quantity of density ρ, for which Q is a source and b is its flux. Over any fixed volume V bounded by the surface ∂V, normal n̂,

\frac{d}{dt} ∫_{V} ρ dV = ∫_{V} Q dV − ∫_{∂V} b · n̂ dS ,

so that the rate of increase of the physical quantity in V is equal to the total rate of generation within V less the amount that flows out of V across the boundary.

Lax-Wendroff schemes for simultaneous conservation equations.

More generally, let u and f[u] be p-vectors satisfying

u_{t} + (f[u])_{x} = 0  (11.2)

Equation (11.2) represents the conservation of p physical quantities. Ideally, we would like our FDM to conserve them also. For example, we could write

U_{n}^{j+1} − U_{n}^{j} = −\frac{1}{2}sΔF_{n}^{j}  where F_{n}^{j} = f [U_{n}^{j}]  and s = k/h .  (11.3)

This scheme is conservative because if we sum over all values of n the right-hand-side all cancels (apart possibly from boundary terms if we are on a finite region of x), so that ∑_{n} U_{n}^{j+1} = ∑_{n} U_{n}^{j}. However, we have seen that (11.3) may be unstable. Alternatively we might rewrite (11.2) in the form

(u_{l})_{t} + ∑_{m=1}^{p} A_{lm}(u_{m})_{x} = 0  for l = 1 … p  where A_{lm} = \frac{∂(f_{l})}{∂(u_{m})} ,  (11.4)

and u_{m} is the m-th component of u. If A[u] is the matrix whose (l, m)th element is A_{lm}, we could then approximate

U_{n}^{j+1} − U_{n}^{j} = −\frac{1}{2}sA (ΔU_{n}^{j})  (11.5)

but this would not be conservative unless A is constant. The Lax-Wendroff idea applied to (11.2) gives

u_{tt} = −(f[u])_{xt} = −(Au_{t})_{x} = (A(f)_{x})_{x}

so that

\begin{aligned}(u_{tt})_{n}^{j} & = \frac{1}{h} [A_{n+1/2}^{j} (f_{x})_{n+1/2}^{j} − A_{n−1/2}^{j} (f_{x})_{n−1/2}^{j}] + O(h^{2}) \\ & = \frac{1}{h^{2}} [A_{n+1/2}^{j} (f_{n+1}^{j} − f_{n}^{j}) − A_{n−1/2}^{j} (f_{n}^{j} − f_{n−1}^{j})] + O(h^{2}) \end{aligned} .  (11.6)

In the above, to the same order of accuracy we can approximate

A_{n+1/2}^{j} ≡ A[u_{n+1/2}^{j}] = A [\frac{1}{2}(u_{n+1}^{j} + u_{n}^{j})] + O(h^{2}) .

Using the estimate (11.6) for u_{tt}, we obtain

U_{n}^{j+1} = U_{n}^{j} − \frac{1}{2}sΔF_{n}^{j} + \frac{1}{2}s^{2}[A_{n+1/2}^{j} (F_{n+1}^{j} − F_{n}^{j}) − A_{n−1/2}^{j} (F_{n}^{j} − F_{n−1}^{j})] . (11.7)

However while (11.7) is is a second order scheme, it still non-conservative if A depends on u, and moreover it is necessary to calculate A_{n+1/2}^{j}. A superior two-step version of the Lax-Wendroff scheme is due to Richtmyer. The system

\begin{matrix}U_{n+1/2}^{j+1/2} & = \frac{1}{2} (U_{n}^{j} + U_{n+1}^{j}) − \frac{1}{2}s (F_{n+1}^{j} − F_{n}^{j}) \\ U_{n}^{j+1} & = U_{n}^{j} − s (F_{n+1/2}^{j+1/2} − F_{n−1/2}^{j+1/2}) \end{matrix} \}  (11.8)

has truncation error O(k^{2}+h^{2}), is conservative and stable if the Courant condition holds. When A is constant, so that F = AU, it reduces to (11.7). For the one-dimensional case p = 1, with f = cu, we obtain the two-step scheme (10.7D) with a different step-size.

The two-step Lax-Wendroff scheme (11.8) may be extended to two space dimensions without difficulty. Consider the equation for u(x, y, t)

u_{t} + f_{x} + g_{y} = 0 and take h_{x} = h_{y} = h .

We define U at the half time-levels for (p, q) = (m, n + 1/2) or (m + 1/2, n) by

U_{pq}^{j+1/2} = \frac{1}{4}[U_{p+,q}^{j} + U_{p−,q}^{j} + U_{p,q+}^{j} + U_{p,q−}^{j}] − \frac{1}{2}s[F_{p+,q}^{j} − F_{p−,q}^{j} + G_{p,q+}^{j} − G_{p,q−}^{j}] ,

where p± denotes p ± 1/2 and similarly q±. At the integer time levels, for (p, q) = (m, n) or (m + 1/2, n + 1/2),

U_{pq}^{j+1} = U_{pq}^{j} − s [F_{p+,q}^{j+1/2} − F_{p−,q}^{j+1/2} + G_{p,q+}^{j+1/2} − G_{p,q−}^{j+1/2}] .  (11.9)

The resulting scheme is second order and conservative.