1.1.3.4. von Neumann Stability Analysis¶

The key to the von Neumann stability analysis is to expand the solution (or error) in a finite Fourier series

1.1.3.4.1. Fourier decomposition of the solution¶

\(\bar{u}_i^n \rightarrow\) exact solution of the discretised equations

\(u_i^n \rightarrow\) numerical solution of the discretised equations (roundoff error, errors in I.C. etc)

\(N \rightarrow\) numerical scheme

\(D \rightarrow\) differential equations

\[u_i^n = \bar{u}_i^n+ \bar{\epsilon}_i^n\]

By definition \(\rightarrow N(\bar{u}_i^n) \equiv 0\)

Apply \(N()\) to equation above

\[N(u_i^n) = N(\bar{u}_i^n + \bar{\epsilon}_i^n)\]

If the Numerical scheme is linear

\[ \begin{align}\begin{aligned}N(u_i^n) = N(\bar{u}_i^n) + N(\bar{\epsilon}_i^n)\\N(u_i^n) = N(\bar{\epsilon}_i^n)\end{aligned}\end{align} \]

If \(N(u_i^n)=0\) represents the numerical scheme
Then the errors satisfy the same equation as the numerical solution, i.e. \(N(\bar{\epsilon}_i^n)=0\)

1.1.3.4.2. Wavenumbers resolvable by a domain¶

Consider a 1D domain \((0, L)\)
Reflect it onto \((-L, 0)\)
Create mesh point, spacing \(\Delta x\)

1.1.3.4.2.1. Maximum wavenumber¶

Shortest resolvable wavelength \(\lambda_{min} = 2 \Delta x\)
Maximum wavenumber \(k_{max} = {{2 \pi} \over {2 \Delta x}} = {\pi \over {\Delta x}}\)

1.1.3.4.2.2. Minimum wavenumber¶

Largest resolvable wavelength \(\lambda_{max} = 2 L\)
Minimum wavenumber \(k_{min} = {{2 \pi} \over {2 L}} = {\pi \over {L}}\)

1.1.3.4.3. Harmonics¶

Mesh index \(i = 0 .... N\), \(x_i = i \Delta x\), \(\Delta x = {L \over N}\)
All harmonics represented in this finite mesh are:

\[k_j = j k_{min} = j {\pi \over L} = j {\pi \over {N \Delta x}}\]

Where: \(j = 0 .... N\) (and \(j = 0\) for a constant solution)

1.1.3.4.4. Phase angle¶

Phase angle is given by:

\[\phi = k_j \Delta x = j {\pi \over N}\]

Covers the whole domain \((-\pi, \pi)\) in steps \(\pi \over N\)

1.1.3.4.5. Decomposition in finite Fourier series¶

\[u_i^n = \sum_{j=-N}^N V_j^n e ^{I k_j x_i} = \sum_{j=-N}^N V_j^n e ^{I k_j i \Delta x} = \sum_{j=-N}^N V_j^n e ^{{I i j \pi} / N}\]

Where \(I = \sqrt{-1}\)
And \(V_j^n\) = amplitude of \(j^{th}\) harmonic

1.1.3.4.6. von Neumann Stability Condition¶

The amplitude of any harmonic many not grow indefinitely in time (as \(n \rightarrow \infty\))
Amplification factor - This is a function of the scheme parameters and \(\phi\) (not n):

\[G = {V^{n+1} \over V^n}\]

Stability condition:

\[ \begin{align}\begin{aligned}\left\vert G \right\vert \le 1\\\forall \phi_j = j {\pi \over N}\\j = -N, N\end{aligned}\end{align} \]

1.1.3.4.7. Example 1: Explicit FD in t, CD in x for 1D linear convection¶

FD in t, CD in x:

\[{{u_i^{n+1} - u_i^n} \over {\Delta t}} + {c \over {2 \Delta x}} ({u_{i+1}^n - u_{i-1}^n})=0\]

Transpose:

\[u_i^{n+1}= u_i^n - {\sigma \over 2} ({u_{i+1}^n - u_{i-1}^n})\]

Replace all terms of the form \(u_{i+m}^{n+k}\) by \(V^{n+k}e^{I(i+m)\phi}\)

\[V^{n+1}e^{Ii\phi} = V^ne^{Ii\phi} - {\sigma \over 2} \left(V^ne^{I(i+1)\phi} - V^ne^{I(i-1)\phi} \right)\]

Divide through by \(e^{Ii\phi}\)

\[V^{n+1} = V^n \left(1 - {\sigma \over 2} \left(e^{I\phi} - e^{-I\phi} \right) \right)\]

Amplification factor - exp function is periodic with period \(2\pi I\), i.e. \(e^{a+bI}=e^a(cosb + Isinb)\):

\[G = {{V^{n+1}} \over {V^n}} = 1-{\sigma \over 2} \left(2I sin \phi \right) = 1- \sigma I sin \phi\]

Stability requires that the “norm of G” be less than 1, i.e. “G” times “G conjugate” is less than 1:

\[GG^* = 1+\sigma^2sin^2\phi \gt 1\]

Hence FD in t, CD in x for 1D linear convection is unconditionally unstable

1.1.3.4.8. Example 2 - Implicit BD in t, CD in x for 1D linear convection¶

BD in t, CD in x:

\[u_i^{n+1}= u_i^n - {\sigma \over 2} ({u_{i+1}^{n+1} - u_{i-1}^{n+1}})\]

\[V^{n+1}= V^n - V^{n+1} {\sigma \over 2} \left({e^{I\phi} - e^{-I\phi}}\right)\]

Amplification factor:

\[G = 1- {\sigma \over 2}G\left({e^{I\phi} - e^{-I\phi}}\right)\]

Re-arrange:

\[G = { 1 \over {1+I\sigma sin \phi}}\]

\[\left\vert G \right\vert^2 = G \cdot G^* = { 1 \over {1+\sigma^2 sin^2 \phi}} \le 1 \qquad \forall \phi\]

BD in t and CD in x for 1D linear convection is unconditionally stable

1.1.3.4.9. Example 3 - Explicit FD in t, BD in x (upwind) for 1D linear convection¶

FD in t, BD in x:

\[u_i^{n+1}= u_i^n - \sigma ({u_{i}^{n} - u_{i-1}^{n}})\]

\[G = {V^{n+1} \over V^n} = 1- \sigma \left(1 - e^{-I\phi} \right) = 1- \sigma + \sigma cos \phi - I \sigma sin\phi\]

Since \(1-cos A=2 sin^2(A/2)\)

\[G = 1- 2\sigma sin^2 \left({\phi \over 2}\right)-I \sigma sin \phi\]

Separate real and imaginary parts of \(G\), \(\xi\), \(\eta\)

\[\xi = Re(G) = 1- 2\sigma sin^2 \left({\phi \over 2}\right) = (1-\sigma)+\sigma cos \phi\]

\[\eta = Im(G) = -\sigma sin \phi\]

Parametric equations for G on complex plane (with \(\phi\) as a parameter) \(\rightarrow\) circle centred on the real axis at point \((1-\sigma)\)

On the complex plane, the stability condition is \(\left\vert G \right\vert \lt 1\)

The circle for \(G\) should be inside the unit circle

Stable for \(0 \le \sigma \le 1\)

Hence FD in t, BD in x for 1D linear convection is conditionally stable

1.1.3.4.10. Example 4 - Implicit BD in t, BD in x for 1D linear convection¶

BD in t, BD in x:

\[u_i^{n+1}= u_i^n - \sigma ({u_{i}^{n+1} - u_{i-1}^{n+1}})\]

\[G = {1 \over {1+\sigma(1-e^{-I\phi}})}\]

Stability:

\[GG^* = {1 \over {(1-\sigma + \sigma cos \phi)^2+\sigma^2sin^2\phi}} \lt 1 \qquad \forall \phi\]

Hence BD in t, BD in x is unconditionally stable

Pattern:

Most explicit schemes are either:

Unconditionally unstable
Conditionally stable

Most implicit schemes are:

Unconditionally stable

Explicit schemes are useless for practical work. However, implicit schemes require more work to solve.

1.1.3.4.11. Example 5 - Explicit FD in t, CD in x for 1D diffusion¶

\[{\partial u \over \partial t} = \nu {\partial^2 u \over \partial x^2}\]

FD in t, CD in x:

\[u_i^{n+1} = u_i^n + \nu {\Delta t \over \Delta x^2}(u_{i+1}^n -2u_i^n+ u_{i-1}^n)\]

\[u_i^{n+1} = u_i^n + \beta(u_{i+1}^n -2u_i^n+ u_{i-1}^n)\]

Amplification factor:

\[G = 1-4 \beta sin^2 \left(\phi \over 2 \right)\]

Stability condition: \(\left\vert 1-4 \beta sin^2 \left(\phi \over 2 \right) \right\vert \le 1\)

Satisfied for:

\[-1 \le 1-4 \beta sin^2 \left(\phi \over 2 \right) \le 1\]

Or:

\[0 \le \beta \le {1 \over 2}\]

Hence BD in t, BD in x is conditionally stable for:

\(\nu \gt 0\) (stability of the physical problem)

and

\({\nu {\Delta t \over \Delta x^2}} \le {1 \over 2}\) (conditional stability of the scheme)

1.1.3.4.12. Example 6 - Implicit BD in t, CD in x for 1D diffusion¶

\[{\partial u \over \partial t} = \nu {\partial^2 u \over \partial x^2}\]

BD in t, CD in x:

\[u_i^{n+1} = u_i^n + \nu {\Delta t \over \Delta x^2}(u_{i+1}^{n+1} -2u_i^{n+1}+ u_{i-1}^{n+1})\]

Transpose:

\[u_i^{n+1} = u_i^n + \beta(u_{i+1}^{n+1} -2u_i^{n+1}+ u_{i-1}^{n+1})\]

Replace all terms of the form \(u_{i+m}^{n+k}\) by \(V^{n+k}e^{I(i+m)\phi}\)

\[V^{n+1}e^{Ii\phi} = V^ne^{Ii\phi} + {\beta} \left(V^{n+1}e^{I(i+1)\phi}- 2V^{n+1}e^{Ii\phi}+ V^{n+1}e^{I(i-1)\phi} \right)\]

Divide through by \(e^{Ii\phi}\)

\[V^{n+1} = V^n + \beta V^{n+1} (e^{I\phi} -2 + e^{-I\phi})\]

Amplification factor - exp function is periodic with period \(2\pi I\), i.e. \(e^{a+bI}=e^a(cosb + Isinb)\):

\[G = {{V^{n+1}} \over {V^n}} = {1 \over {1-2\beta(cos\phi -1)} }\]

Stability requires that the magnitude of G is less than 1:

\[ \begin{align}\begin{aligned}cos \phi = 1 \qquad \left\vert G \right\vert = \left\vert {1 \over {1-2\beta} } \right\vert\\cos \phi = 0 \qquad \left\vert G \right\vert = \left\vert {1 \over {1+2\beta} } \right\vert\\cos \phi = -1 \qquad \left\vert G \right\vert = \left\vert {1 \over {1-4\beta} } \right\vert\\\beta \rightarrow \infty \qquad \left\vert G \right\vert \rightarrow 0\\\beta = 0 \qquad \left\vert G \right\vert = 1\\\beta \rightarrow -\infty \qquad \left\vert G \right\vert \rightarrow 0\end{aligned}\end{align} \]

\[\left\vert G \right\vert \le 1\]

Hence BD in t, CD in x for 1D diffusion is unconditionally stable

1.1.3.4.13. Summary¶

Stability conditions place a limit on the time step for a given spatial step.
This has a physical interpretation - the solution progresses too rapidly in time - especially a problem for convection dominated flows and compressible flows at the speed of sound - if c is large \(\Delta t\) must be small.
Implicit schemes have no limit on timestep. Implicit vs explicit is a debatable area for different applications.
For the diffusion equation, the explicit time step restriction here is not too severe. But numerical diffusion can be large, depending on \(\Delta x\).
The stability of linear schemes is well understood. But we have also neglected boundaries.