Computational Fluid Dynamics - Elementary (h)

Stability of the Finite Difference Schemes

Von Neumann method

A very versatile tool for analysing stability is the Fourier method developed by von Neumann. In this method, the initial values at mesh points are expressed in terms of finite Fourier series, and we consider the growth of individual Fourier components. In this method, it is not required to find eigenvalues or matrix norms. Therefore, it is one of the simplest ways to evaluate stability.

A finite sine or cosine series expansion in the interval $ x \in {a, b} $ takes the form:
\begin{equation} \sum\limits_s a_s \sin \frac{s \pi x}{L} \ \ \ , \ \ \ \sum\limits_s b_s \cos \frac{s \pi x}{L} \end{equation} where $ L = b - a $.

Now consider an individual component written in complex exponential form at a mesh point $ x = x_j = a +j\Delta x $:
\begin{equation} A_s e^{\frac{isx\pi}{L}} = A_s e^{\frac{isa\pi}{L}} e^{\frac{isj \Delta x\pi}{L}} = \bar{A}_s e^{i\alpha_sj \Delta x}; \end{equation} where $ \alpha _s = s \pi / L $.

Given initial data we can express the initial values as:
\begin{equation} w_j^0 = \sum\limits_{S=0}^n \bar{A}_s e^{i\alpha_sj \Delta x} \ \ \ \ \ \ J = 0, 1, \dots, n. \end{equation} We use the $ n+1 $ equations above to determine the $ n+1 $unknowns $ \bar{A} $. Now we want to find out how each Fourier mode develops in time. So we assume a simple separable solution of the form:
\begin{equation} w_j^k = \sum\limits_{S=0}^n \bar{A}_s e^{i\alpha_sj \Delta x} e^{\Omega t_k} = \sum\limits_{S=0}^n \bar{A}_s e^{i\alpha_sj \Delta x} e^{\Omega k \Delta t} = \sum\limits_{S=0}^n \bar{A}_s e^{i\alpha_sj \Delta x} \xi^k; \end{equation} where $ \xi = e^{\Omega \Delta t} $ which is called the amplification factor. For stability we thus require:
\begin{equation} | \xi | \leq 1. \end{equation}

Example

As a simple example, consider the 1-Dimensional heat equation:
\begin{equation} \frac{\partial u}{\partial t} = k \frac{\partial^2 u}{\partial x^2}. \end{equation} This differential equation in discretized form can be written as:
\begin{equation} u_j^{n+1} = u_j^n + r(u_{j+1}^n - 2u_j^n + u_{j-1}^n); \end{equation} where $ r = \frac{k \Delta t}{\Delta x^2} $. For each term we can write the Fourier mode developed in time:
\begin{equation} \begin{aligned} & w_j^n = e^{i \alpha_s j \Delta x}\xi^n; \\ & w_j^{n+1} = e^{i \alpha_s j \Delta x}\xi^{n+1}; \\ & w_{j+1}^n = e^{i \alpha_s (j+1) \Delta x}\xi^n; \\ & w_{j-1}^n = e^{i \alpha_s (j-1) \Delta x}\xi^n. \end{aligned} \end{equation} These terms should also satisfy the heat equation:
\begin{equation} w_j^{n+1} = w_j^n + r(w_{j+1}^n - 2w_j^n + w_{j-1}^n). \end{equation} Now we can plug in our defined single separable solutions into the discretized equation as follows (after simplification):
\begin{equation} \xi = 1 + r(e^{i \alpha_s \Delta x} + e^{-i \alpha_s \Delta x} - 2). \end{equation} Using the following definitions:
\begin{equation} \cos(\alpha_s \Delta x) = \frac{e^{i \alpha_s \Delta x} + e^{-i \alpha_s \Delta x}}{2} \ \ \ , \ \ \ \sin^2 \left(\frac{\alpha_s \Delta x}{2}\right) = \frac{1 - \cos(\alpha_s \Delta x)}{2}; \end{equation} the above equation can be written as:
\begin{equation} \xi = 1 - 4r \sin^2 \left(\frac{\alpha_s \Delta x}{2}\right). \end{equation} The necessary and sufficient condition for the error to remain bounded is that $ | \xi | \leq 1 $, so:
\begin{equation} \left|1 - 4r \sin^2 \left(\frac{\alpha_s \Delta x}{2}\right)\right| \leq 1; \end{equation} which has to be hold for all values of $ sin^2(\dots) $, therefore:
\begin{equation} \frac{k \Delta t}{\Delta x^2} \leq \frac{1}{2}; \end{equation} represents the stability requirement of this finite difference scheme.