Observation

When using the upwind scheme for the solution of advection equations, there is a critical timestep size, $\displaystyle \Delta t$ , above which the solution becomes unstable.

Problem Statement

Use von-Neumann’s stability analysis to establish the timestep-size

requirements for the upwind scheme to generate stable solutions.

Approach

To solve this problem, we need to do the following:

Consider a 1-D advection equation with a constant velocity.
Discretize the 1-D advection equation using forward-Euler in time and upwind in space.
Given the difference equation, introduce an artificial perturbation as an initial condition.
Estimate the growth/decay properties of artificial perturbation and derive the stability requirement for this problem.

Show Me the Math

Consider the advection of a quantity $\displaystyle \phi$ at constant speed $\displaystyle v$ .

This motion can be described by the advection equation

$\displaystyle \frac{\partial\phi}{\partial t}=-v\frac{\partial\phi}{\partial x}$

We assume that $\displaystyle v>0$ for simplicity of the exhibition. Our purpose is to study the stability of a forward-Euler, upwind discretization of this equation. For this purpose, we set

$\displaystyle \frac{\phi_{i}^{n+1}-\phi_{i}^{n}}{\Delta t}=-v\frac{\phi_{i}^{n}-\phi_{i-1}^{n}}{\Delta x}$

$\displaystyle \phi_{i}^{n+1}=\phi_{i}^{n}-v\frac{\Delta t}{\Delta x}(\phi_{i}^{n}-\phi_{i-1}^{n})$

Here, I will assume that you are familiar with von-Neumann stability analysis. If not, then see my description \ref{sub:Von-Neumann-Stability}. Without further ado, we substitute

$\displaystyle \phi_{i}^{n}\equiv\sum_{m}e^{\alpha t}e^{ik_{m}x}$

We then have

$\displaystyle \sum_{m}e^{\alpha(t+\Delta t)}e^{ik_{m}x}=\sum_{m}e^{\alpha t}e^{ik_{m}x}-v\frac{\Delta t}{\Delta x}\left[\sum_{m}e^{\alpha t}e^{ik_{m}x}-\sum_{m}e^{\alpha t}e^{ik_{m}(x-\Delta x)}\right]$

but, by virtue of linearity, we can safely drop the summation and

investigate the behaviour of a single mode. In other words

$\displaystyle e^{\alpha(t+\Delta t)}e^{ik_{m}x}=e^{\alpha t}e^{ik_{m}x}-v\frac{\Delta t}{\Delta x}\left[e^{\alpha t}e^{ik_{m}x}-e^{\alpha t}e^{ik_{m}(x-\Delta x)}\right]$

then, dividing by $\displaystyle e^{\alpha t}e^{ik_{m}x}$ , we have

$\displaystyle A\equiv e^{\alpha\Delta t}=1-C\left[1-e^{-ik_{m}\Delta x}\right];\quad C\equiv v\frac{\Delta t}{\Delta x}$

where $\displaystyle A$ is the error amplification factor. Using the exponential

trigonometric identities, we have

$\displaystyle A=1-C\left[1-\cos(k_{m}\Delta x)\right]-iC\sin(k_{m}\Delta x)$

The method is stable is the amplification factor is less than 1. The

magnitude of the amplification factor is

$\displaystyle A^{2}=\left\{ 1-C\left[1-\cos(k_{m}\Delta x)\right]\right\} ^{2}+C^{2}\sin^{2}(k_{m}\Delta x)$

$\displaystyle A^{2}=1-2C\left[1-\cos(k_{m}\Delta x)\right]+C^{2}\left[1-\cos(k_{m}\Delta x)\right]^{2}+C^{2}\sin^{2}(k_{m}\Delta x)$

but

$C^{2}\left[1-\cos(k_{m}\Delta x)\right]^{2}+C^{2}\sin^{2}(k_{m}\Delta x) =C^{2}\left[1-2\cos(k_{m}\Delta x)+\cos^{2}(k_{m}\Delta x)\right]+C^{2}\sin^{2}(k_{m}\Delta x)$