A Simple Algorithm

As a first attempt to solve (2.9) we consider using centred differences for the space derivative and Euler's method for the time part

$$u_j^{n+1} = u_j^n - {c \delta t\over 2\delta x} \left(u_{j+1}^n - u_{j-1}^n\right)$$ (2.11)

where the subscripts ${}_j$ represent space steps and the superscripts ${}^n$ time steps. By analogy with the discussion of the Euler and Leap-Frog methods we can see immediately that this method is 1st order accurate in $t$ and 2nd order in $x$. We note firstly that

$$u_j^{n+1} - u_j^n \approx \delta t \left.\partial u\over\... ... \left.\partial^2 u\over\partial t^2\right\vert _j^n + \ldots$$ (2.12)

whereas

$$u_{j+1}^n - u_{j-1}^n \approx 2\delta x \left.\partial u\... ... \left.\partial^3 u\over\partial x^3\right\vert _j^n + \ldots$$ (2.13)

and substitute these forms into (2.11) to obtain

$$\delta t \left.\partial u\over\partial t\right\vert _j^n +... ...\partial^3 u\over\partial x^3\right\vert _j^n + \ldots\right.$$ (2.14)

so that, when the original differential equation (2.9) is subtracted, we are left with a truncation error which is 2nd order in the time but 3rd order in the spatial part. The stability is a property of the time, $t$, integration rather than the space, $x$. We analyse this by considering a plane wave solution for the $x$-dependence by substituting $u_j^n = v^n\exp(\i k x_j)$ to obtain

$$v^{n+1} e^{\i k x_j} = v^n e^{\i k x_j} - {c \delta t\over 2\delta x} v^n \left(e^{i k x_{j+1}} - e^{\i k x_{j-1}}\right)$$ (2.15)

or, after dividing out the common exponential factor,

$$v^{n+1} = \left[ 1 - \i{c\delta t\over\delta x} \sin(k \delta x) \right]v^n.$$ (2.16)

Since the wave and advection equations express a conservation law the solution should neither grow nor decay as a function of time. If we substitute $v^n \mapsto v^n + \delta v^n$ and subtract (2.16) we obtain an equation for $\delta v^n$

$$\delta v^{n+1} = \left[ 1 - \i{c\delta t\over\delta x} \sin(k \delta x) \right]\delta v^n$$ (2.17)

which is identical to (2.16). Hence the stability condition is simply given by the quantity in square brackets. We call this the amplification factor and write it as

$$g = 1 - \i\alpha$$ (2.18)

where

$$\alpha = {c\delta t\over\delta x}\sin(k \delta x).$$ (2.19)

As $g$ is complex the stability condition becomes

$$\left\vert g\right\vert^2 = 1 + \alpha^2 \le 1 \mbox{{} for all $k$\ {}}.$$ (2.20)

The condition must be fulfilled for all wave vectors $k$; otherwise a component with a particular $k$ will tend to grow at the expense of the others. This is known as the von Neumann stability condition^2.2. Unfortunately it is never fulfilled for the simple method applied to the advection equation: i.e. the method is unstable for the advection equation.