4. Suggested solutions weekly assignments - lecture 13#

Only answers to the stability computations are shown below.

Find the absolute stability limit of the diffusion equation discretized with a Crank-Nicolson scheme. Use both a finite difference method \( \boldsymbol{u}^{n+1} - \boldsymbol{u}^n = \frac{\Delta t}{2}D^{(2)}(\boldsymbol{u}^{n+1}+\boldsymbol{u}^n) \) and a comparable Galerkin method.

Insert for \(\boldsymbol{u}^{n}=g^n \boldsymbol{u}^0\) to get

\[ g^n(g-1) \boldsymbol{u}^0 = \frac{\Delta t}{2}D^{(2)}g^n(g+1) \]

Use the eigenvalues of the differentiation matrix: \(D^{(2)}\boldsymbol{u}^0 = \lambda \boldsymbol{u}^0\)

\[ g^n(g-1) \boldsymbol{u}^0 = \frac{\Delta t}{2}\lambda \boldsymbol{u}^0 g^n(g+1) \]

Divide by \(g^n \boldsymbol{u}^0\) to get

\[ g-1 = \frac{\Delta t}{2} \lambda (g+1) \]

Rearrange to get

\[ g = \frac{2+\Delta t \lambda}{2-\Delta t \lambda} \]

Compute the eigenvalues:

import numpy as np
from scipy import sparse

N = 11
D2 = sparse.diags((1, 1, -2, 1, 1), (-N, -1, 0, 1, N), shape=(N+1, N+1))
Lambda = np.linalg.eig(D2.toarray())[0]
print(Lambda)

[-4.00000000e+00 -3.73205081e+00 -3.00000000e+00 -2.00000000e+00
 -1.00000000e+00 -4.53075522e-16 -2.67949192e-01 -2.67949192e-01
 -1.00000000e+00 -3.73205081e+00 -2.00000000e+00 -3.00000000e+00]

All eigenvalues are real and negative. The smallest value is \(\lambda = -4\). Insert into equation for \(g\)

\[ g = \frac{2-4 \Delta t }{2+4\Delta t} \]

Obviously, \(|g|\le 1\) for all \(\Delta t \ge 0\), so the Crank-Nicolson method is unconditionally stable for the diffusion equation.

Consider the leap-frog finite difference scheme \( \boldsymbol{u}^{n+1} - \boldsymbol{u}^{n-1} = {2\Delta t}D^{(2)}\boldsymbol{u}^n, \) suggested by Richardson in 1910. Explain why this scheme is unconditionally unstable. To this end you can use that all the eigenvalues of \(D^{(2)}\) are real.

Insert for \(\boldsymbol{u}^n = g^n \boldsymbol{u}^0\) to get

\[ g^n(g - g^{-1}) \boldsymbol{u}^0 = 2 \Delta t D^{(2)} g^n \boldsymbol{u}^0 \]

Use eigenvalues and divide by \(g^n \boldsymbol{u}^0\) to get

\[ g - g^{-1} = 2 \Delta t \lambda \]

Solve quadratic equation for \(g\)

\[ g = \Delta t \lambda \pm \hat{\imath} \sqrt{1+(\Delta t \lambda)^2} \]

And since all \(\lambda\) are real we find that \(|g|=\sqrt{1+2(\Delta t \lambda)^2} \ge 1 \, \forall \, \Delta t \ge 0\). Hence the leap-frog scheme is unconditionally unstable!