The one-dimensional advection equation and Upwinding

Consider the problem for u(x, t)

u_{t} + c u_{x} = 0  for t > 0,  in u(x, 0) = f(x) ,  (10.5)

where c is a positive constant. The characteristics of this equation are dx/dt = c or x − ct = constant. Furthermore, u is constant along each characteristic and so the solution is u = f(x − ct). This represents a wave travelling in the positive x-direction with speed c, without change of size or shape. Can we develop FDMs with these conservation properties? For future reference, we observe that for a Fourier mode with f(x) = exp(ixξ/h)

u_{n}^{j+1} ≡ u(nh, (j + 1)k) = u_{n}^{j} exp(−iξq)  where q = ck/h > 0  (10.6)

is the Courant number. The real solution therefore has a growth factor λ = exp(−iξq). We note that |λ| = 1 so that there is no growth or decay. We can model (10.5) in many different ways. We shall use a forward difference for u_{t} and consider 6 schemes:

\frac{U_{n}^{j+1} − U_{n}^{j}}{−q} = \{ \begin{matrix}A & \frac{1}{2}ΔU_{n}^{j} ≡ \frac{1}{2}(U_{n+1}^{j} − U_{n−1}^{j}) & Explicit, centred \\ B & U_{n+1}^{j} − U_{n}^{j} & Explicit, forwards \\ C & U_{n}^{j} − U_{n−1}^{j} & Explicit, backwards \\ D & \frac{1}{2}ΔU_{n}^{j+1/2} & Two-step, centred \\ E & \frac{1}{2} [\frac{1}{2}ΔU_{n}^{j+1} + \frac{1}{2}ΔU_{n}^{j}] & Crank-Nicolson \\ F & \frac{1}{2}ΔU_{n}^{j} − \frac{1}{2}qδ^{2}U_{n}^{j} & Lax-Wendroff \end{matrix}  (10.7)

(See also the ‘Box scheme,’ on example sheet 2.) We use the Fourier method to investigate the stability of all of these schemes, looking for a solution U_{n}^{j} = λ^{j} exp(inξ). We find:

λ = \{ \begin{matrix}A & 1 − iq sin ξ \\ B & 1 − q(e^{iξ} − 1) \\ C & 1 − q(1 − e^{−iξ}) \\ D & \frac{1}{4} [iq sin ξ ± \sqrt{4 − q^{2} sin^{2} ξ}]^{2}  (two steps) \\ E & \frac{2 − iq sin ξ}{2 + iq sin ξ} \\ F & 1 − iq sin ξ − 2q^{2} sin^{2} \frac{1}{2}ξ \end{matrix} .  (10.8)

We recall that if |λ| > 1 + O(k) the method is unstable; if |λ| < 1 it is dissipative, while if |λ| = 1 is it conservative. We see therefore that

|λ|^{2} = \{ \begin{matrix}A & 1 + q^{2} sin^{2} ξ & stable only if q^{2} = O(k) i.e. k ∼ h^{2} \\ B & 1 + 2q(1 + q)(1 − cos ξ) & ⩾ 1 ∀ q, ξ  hopelessly unstable \\ C & 1 − 2q(1 − q)(1 − cos ξ) & ⩽ 1 ∀ q ⩽ 1  stable, dissipative \\ D & 1 & provided q^{2} sin^{2} ξ ⩽ 4  or \frac{1}{2}q ⩽ 1 \\ E & 1 & ∀ q  stable and conservative \\ F & 1 − 4q^{2}(1 − q^{2}) sin^{4} \frac{1}{2}ξ & ⩽ 1 if q ⩽ 1 dissipative, but less than C \end{matrix}

We may also be interested in arg λ, which determines the phase of each mode. In case E, for example,

arg λ = −2 tan^{−1}(\frac{1}{2}qξ) ≈ −qξ when ξ is small.

All the above schemes agree well with the exact solution, for which arg λ = −qξ, when ξ is small. The waves with higher values of ξ are dispersive; their phase velocity varies with frequency. The computer demonstration will illustrate the effects of this. Although the amplitude of each wave component may be conserved, phase differences develop which alter the shape of the whole.

We can see from the above the importance of the CFL condition, q ⩽ 1. Another very important conclusion we can draw is that backwards, or UPWIND differences are a good idea. Case C is stable provided the CFL condition holds, whereas B is unconditionally unstable. If c < 0, the stable scheme is B. In general, the Upwind scheme can be written

U_{n}^{j+1} = U_{n}^{j} − scΔU_{n}^{j} + s|c|δ^{2}U_{n}^{j}  where s = \frac{k}{2h} .  (10.9)

The need for upwind differences can be interpreted in terms of the characteristics of (10.5), x − ct = constant which must pass through the Numerical Domain of Dependence.

Upwinding for simultaneous equations

Let u(x, t) be a p-dimensional vector satisfying the equation u_{t} + Au_{x} = d where A is a p × p matrix, which for simplicity we shall assume to be constant. We shall investigate how the upwinding method generalises to this problem, by diagonalising the matrix A, which we assume to have p distinct eigenvalues λ_{1} … λ_{p} and corresponding eigenvectors. We make the linear transformation u = Sv where S is the matrix whose columns are the eigenvectors of A, so that

S^{−1}AS = D ≡ diag(λ_{1}, λ_{2}, …, λ_{p}) and so A = SDS^{−1} .

Then v_{i}, the i-th component of v, satisfies the equation

(v_{i})_{t} + λ_{i}(v_{i})_{x} = (S^{−1}d)_{i} .

We have thus separated the problem into p separate ones, each with its own characteristic family, dx/dt = λ_{i}, and we can use upwinding on each one in turn. Note that if any of the λ_{i} are complex, then some of the p equations are elliptic, and we cannot use a time-stepping approach for them. We assume here that all the λ_{i} are real. From (10.9), the upwind scheme for v is

V_{n}^{j+1} = V_{n}^{j} − sDΔV_{n}^{j} + sD^{+}δ^{2}V_{n}^{j} + kS^{−1}d_{n}^{j} ,

where D^{+} ≡ diag(|λ_{1}|, |λ_{2}|, …, |λ_{p}|), and V is the FDA to v. Then transforming back, defining U = SV so that V = S^{−1}U, we find

U_{n}^{j+1} = U_{n}^{j} − sAΔU_{n}^{j} + sA^{+}δ^{2}U_{n}^{j} + kd_{n}^{j}  where A^{+} = SD^{+}S^{−1} .  (10.10)

This is the required generalisation of upwinding for p simultaneous equations. If all the eigenvalues of A are positive, then D^{+} = D and A^{+} = A, while A^{+} = −A if they are all negative. Otherwise, A^{+} bears no simple relation to A, and we must be very careful how we define “up” when upwinding. The stability condition for (10.10) is

\frac{k}{h} max_{i=1...p} |λ_{i}| ⩽ 1 .  (10.11)