Theory exercise 5

Fredrik Eikeland Fossan

Mar 18, 2019

Problem 1: Stability

In this problem we will consider the diffusion equation $$ \begin{align} \frac{\partial u}{\partial t} = \alpha\frac{{\partial}^2u}{\partial x^2} \,, \label{eq:diffusion} \end{align} $$

and the advection equation: $$ \begin{align} \frac{\partial u}{\partial t} + a_0\frac{\partial u}{\partial x} = 0 \,, \quad a_0 > 0 \,. \label{eq:advection} \end{align} $$

Notation: $$ \begin{align} t_n = \Delta t \cdot n \,, n=0,1\,\dots N \,, \quad x_j = \Delta x \cdot j \,, j=0,1\,\dots K \,. \label{eq:notation} \end{align} $$

a) Discretize Eq. \eqref{eq:diffusion} with first order backward differences for $ \frac{\partial u}{\partial t} $ and central differences for $ \frac{{\partial}^2u}{\partial x^2} $ and show that resulting discretized equation takes the form: $$ \begin{align} u_j^{n + 1} = u_j^{n} + D\left(u_{j + 1}^{n + 1} - 2u_{j}^{n+1} + u_{j - 1}^{n + 1}\right) \,, \label{eq:Laasonen} \end{align} $$

and write out the expression for $ D $. Is the scheme explicit or implicit?

Answer.

$ D= \alpha \frac{\Delta t}{{\Delta x}^2} $. The scheme is implicitt since a solution of a system of algebraic equations is necessary to obtain the unknowns of the next time-step.

b) Use von Neumann stability analysis and show that Eq. \eqref{eq:Laasonen} is stable for all values of $ D $ (unconditionally stable).

Answer.

We introduce $$ \begin{align*} E_j^n = G^n e^{i \cdot \beta \cdot x_j} = G^n e^{i \cdot \beta \cdot \Delta x \cdot j} = G^n e^{i \cdot \delta \cdot j} \, \quad \delta = \beta \cdot \delta x \,. \end{align*} $$ Substitution of $ E_j^{n + 1}, E_{j+1}^{n + 1}, E_{j-1}^{n + 1}, E_j^n $ into Eq. \eqref{eq:Laasonen} we obtain $$ \begin{align*} G^{n+1} e^{i \cdot \delta \cdot j} = G^n e^{i \cdot \delta \cdot j} + D\left(G^{n+1} e^{i \cdot \delta \cdot \left(j + 1\right)} - 2G^{n+1} e^{i \cdot \delta \cdot j} + G^{n+1} e^{i \cdot \delta \cdot \left(j -1\right)}\right) \,, \end{align*} $$ and further by division of $ G^n e^{i \cdot \delta \cdot j} $ $$ \begin{align*} G = 1 + D\left(G e^{i \cdot\delta} - 2G + G e^{-\,i \cdot\delta}\right) = 1 -2 \, D \cdot G \cdot \left(1 - cos \, \delta\right) = 1 - 4\, G \cdot D \,sin^2 \,\frac{\delta}{2} \rightarrow G = \frac{1}{1 + 4 \,D\,sin^2\,\frac{\delta}{2}}\,, \end{align*} $$ from which it follows that $ \lvert G \rvert \leq 1 $ for any value of $ D $, and the scheme is unconditionally stable.

c) Another discretized version (scheme) for Eq. \eqref{eq:diffusion} is given by: $$ \begin{align} D\,u_{j - 1}^{n + 1} - \left(2\,D+3\right)u_{j}^{n+1} + D\,u_{j + 1}^{n + 1} = -2D\,u_{j - 1}^{n} + \left(4\,D-3\right)u_{j}^{n} - 2D\,u_{j + 1}^{n}\,. \label{eq:other} \end{align} $$

Use von Neumann stability analysis and show that the the amplification factor $ G $, may be expressed as $$ \begin{align} G =\frac{3 - 8\,D\cdot sin^2\frac{\delta}{2}}{3 + 4\,D\cdot sin^2\frac{\delta}{2}}\,, \quad \delta = \beta \cdot \Delta x \,, \label{eq:other_amplification} \end{align} $$ and use this to find an expression (in terms of $ D $) for the stability limit of the scheme.

Answer.

Substitution of $ E_j^{n + 1}, E_{j+1}^{n + 1}, E_{j-1}^{n + 1}, E_j^n $ into Eq. \eqref{eq:other} and by division of $ G^n e^{i \cdot \delta \cdot j} $, we obtain $$ \begin{align*} D\,G e^{-i \cdot \delta} - \left(2\,D+3\right)G + D\,G e^{i \cdot \delta} = -2D\,e^{-i \cdot \delta} + \left(4\,D-3\right) - 2D\,e^{i \cdot \delta}\,,\\ \\ \rightarrow G \left[D\left(e^{i \, \delta} + e^{-i \, \delta}\right) - 2\,D - 3\right] = -2\,D\left(e^{i \, \delta} + e^{-i \, \delta}\right) + 4\,D -3\\ \rightarrow G \left[2\,D\left(cos \, \delta -1) - 3\right)\right] = 4\,D\left(1 - cos \, \delta\right) - 3 \\ \rightarrow G \left[-2\,D\cdot 2 \, sin^2 \frac{\delta}{2} - 3\right] = 4\,D\cdot 2 \, sin^2 \frac{\delta}{2} - 3 \\ \rightarrow G = \frac{3 - 8\,D \cdot sin^2 \frac{\delta}{2}}{3 + 4\,D\cdot sin^2 \frac{\delta}{2}} \end{align*} $$ The strict von Neumann condition gives $$ \begin{align*} \lvert G \rvert \leq 1 \rightarrow -1 \leq \frac{3 - 8\,D \cdot sin^2 \frac{\delta}{2}}{3 + 4\,D\cdot sin^2 \frac{\delta}{2}} \leq 1 \,, \end{align*} $$ since G is real. The right hand side is satisfied for any $ D $. The left hand side gives $$ \begin{align*} \lvert G \rvert \leq 1 \rightarrow -1 \leq \frac{3 - 8\,D \cdot sin^2 \frac{\delta}{2}}{3 + 4\,D\cdot sin^2 \frac{\delta}{2}} \rightarrow D \leq \frac{3}{2\,sin^2 \frac{\delta}{2}} \,, \end{align*} $$ which gives the stability limit $ D \leq \frac{3}{2} $.

d) Forward in time central in space (FTCS) discretization of Eq. \eqref{eq:advection} results in the following scheme: $$ \begin{align} u_j^{n + 1} = u_j^{n} - \frac{C}{2}\left(u_{j + 1}^{n} - u_{j- 1}^{n}\right) \,, \label{eq:FTCS} \end{align} $$ where $ C=a_0 \cdot \frac{\Delta t}{\Delta x} $. The FTCS scheme is unconditionally unstable, but a modification of Eq. \eqref{eq:FTCS}, suggested by Peter Lax, is to replace $ u_j^n $ by $ u_j^n = \frac{1}{2}\left(u_{j+1}^n + u_{j-1}^n\right) $. Introduce this modification and show that the resulting scheme i stable for $ C \leq 1 $, by use of the criterion for positive coefficients (PC-criterion).

Answer.

The resulting scheme may be expressed as $$ \begin{align*} u_j^{n + 1} = \frac{1}{2}\left[\left(u_{j+1}^n + u_{j-1}^n\right) - C \cdot \left(u_{j+1}^n - u_{j-1}^n\right)\right] = \frac{1}{2}\left[\left(1 - C\right)\,u_{j+1}^n + \left(1 + C\right)u_{j-1}^n\right]\,, \end{align*} $$ where we have $ a_1 = \left(1-C\right)/2 $ and $ a_2 = \left(1+C\right)/2 $ which gives $ a_1 + a_2 = 1 $. The criterion for positive coefficients further requires that all coefficients are positive. with $ C > 0 $, $ a_2 $ is always positive. Finally $ a_1 \geq 0 $ requires that $ C \leq 1 $.

e) Use von Neumanns stability analysis on and show that the amplification factor $ G $, for the resulting scheme in the previous sub-exercise, may be expressed as $$ \begin{align} G = cos \, \delta - i \cdot C \cdot sin \,\delta \,, \label{eq:Lax_amplification} \end{align} $$

and show that you obtain the same stability limit as with the PC-criterion.

Answer.

Substitution of $ E_j^{n + 1}, E_{j+1}^{n}, E_{j-1}^{n}, E_j^n $ and by division of $ G^n e^{i \cdot \delta \cdot j} $, we obtain $$ \begin{align*} G = \frac{1}{2}\left[\left(1 - C\right)\,e^{i \cdot \delta} + \left(1 + C\right)e^{-i \cdot \delta}\right] = cos \, \delta - i \cdot C \cdot sin \, \delta\,, \end{align*} $$

which gives $$ \begin{align*} {\lvert G \rvert}^2 = cos^2 \, \delta + C^2 \, sin^2 \, \delta = 1 - sin^2\, \delta + C^2 \, sin^2 \, \delta \leq 1 \rightarrow sin^2 \, \delta \left(1 - C^2\right) \geq 0\,, \end{align*} $$ and requires that $ C \leq 1 $.

f) Discretize Eq. \eqref{eq:advection} with first order forward differences for both the temporal and spatial derivatives, and use von Neumann stability analysis to show that the resulting scheme is unstable for all values of C.

Hint.

Show that the amplification factor may be written as: $$ \begin{align} G = 1 + C \cdot \left(1 - cos \, \delta \right) - C \cdot i \cdot sin \, \delta \,, \quad \rightarrow \quad \lvert G \rvert^2 = 1 + 2 \cdot C \cdot \left(1 - cos \, \delta\right) \cdot \left(1 + C\right) \,, \label{eq:ftfs_amplification} \end{align} $$ and introduce the substitution $ 1 - cos \, \delta = 2 \cdot sin^2 \frac{\delta}{2} $

Answer.

The resulting scheme may be expressed as $$ \begin{align*} u_j^{n + 1} = u_j^n - C\,\left(u_{j+1}^n - u_j^n\right)\,. \end{align*} $$ Substitution of $ E_j^{n + 1}, E_{j+1}^{n}, E_j^n $ and by division of $ G^n e^{i \cdot \delta \cdot j} $, we obtain $$ \begin{align*} G = 1 - C \cdot\left(e^{i \cdot \delta} - 1\right) = 1 + C \cdot \left(1 - cos \, \delta\right) - i \cdot C \cdot sin \, \delta\,, \end{align*} $$

which gives $$ \begin{align*} {\lvert G \rvert}^2 = \left[1 + C \cdot \left(1 - cos \, \delta \right)\right]^2 + C^2 \, sin^2 \, \delta \\ = 1 + 2C - 2C \, cos \, \delta + 2C^2 - 2C^2 cos \, \delta \\ = 1 + 2C \cdot \left(1 - cos \, \delta \right) + 2C^2 \cdot \left(1 - cos \, \delta \right) \\ = 1 + 2C \cdot \left(1 - cos \, \delta \right) \cdot \left(1 + C\right) \\ = 1 + 4C \cdot sin^2 \, \frac{\delta}{2} \cdot \left(1 + C\right) \end{align*} $$ and we see that $ \lvert G \rvert \geq 1 $ for all values $ C $. $ G = 1 $ for $ C = 0 $, however this requires a time-step of 0.