Problem Set 5 (Numeric PDEs)

1

Consider the 1D wave equation

$\frac{\partial^2 u}{\partial t^2} = v^2 \frac{\partial^2 u}{\partial x^2}$

(a)

Write down the straightforward finite-difference approximation.

$\frac{u_j^{n + 1} - 2u_j^n + u_j^{n - 1}}{(\Delta t)^2} = v^2 \frac{u_{j+1}^n - 2u_j^n + u_{j-1}^n}{(\Delta x)^2}$

(b)

What order approximation is this in time and in space?

Second order in space and time.

(c)

Use the von Neumann stability criterion to find the mode amplitudes.

Start with an ansatz of $u_j^n = A^n e^{ikj \Delta x}$ . (Here the “exponent” on $u$ is an index, while the exponent on $A$ is an exponent.) Note that $e^{ix} + e^{-ix} = 2 \cos(x)$ (which can be proven by looking at the Taylor expansions), and $\cos(x) = -2\sin(x/2)^2 + 1$ . Then

$\begin{aligned} \left(A^{n+1} - 2A^n + A^{n-1} \right) e^{ikj \Delta x} &= \frac{v^2 (\Delta t)^2}{(\Delta x)^2} A^n \left(e^{ik(j + 1) \Delta x} -2e^{ikj \Delta x} + e^{ik(j - 1) \Delta x} \right) \\ A - 2 + \frac{1}{A} &= \frac{v^2 (\Delta t)^2}{(\Delta x)^2} \left(e^{ik \Delta x} - 2 + e^{-ik \Delta x} \right) \\ &= \frac{2 v^2 (\Delta t)^2}{(\Delta x)^2} \left(\cos(k \Delta x) - 1 \right) \\ &= \frac{-4 v^2 (\Delta t)^2}{(\Delta x)^2} \sin \left(\frac{k \Delta x}{2} \right)^2 \\ 0 &= A^2 + 2 \left( \frac{2 v^2 (\Delta t)^2}{(\Delta x)^2} \sin \left(\frac{k \Delta x}{2} \right)^2 - 1 \right) A + 1 \\ &= A^2 + 2 \alpha A + 1 \end{aligned}$

where $\alpha$ is just defined to avoid rewriting that big subexpression. So by the quadratic formula, the two possible amplitudes for mode $k$ are

$\begin{aligned} A_\pm &= \frac{-2 \alpha \pm \sqrt{4 \alpha^2 - 4}}{2} \\ &= -\alpha \pm \sqrt{\alpha^2 - 1} \end{aligned}$

(e)

Use this to find a condition on the velocity, time step, and space step for stability. (Hint: consider the product of the two amplitude solutions.)

$\begin{aligned} A_+ A_- &= \left( -\alpha + \sqrt{\alpha^2 - 1} \right) \left( -\alpha - \sqrt{\alpha^2 - 1} \right) \\ &= \alpha^2 - (\alpha^2 - 1) \\ &= 1 \end{aligned}$

Thus the magnitudes of the two amplitudes (we never said they were real) have to multiply to one. This means that one will be unstable unless both are exactly equal to one.

If $\alpha^2 \geq 1$ , then both amplitudes are real. To be stable, both must be equal to one, which can only be true if $\sqrt{\alpha^2 - 1} = 0$ , i.e. $\alpha = \pm 1$ .

If $\alpha^2 < 1$ , then the amplitudes are $A_\pm = -\alpha \pm i \sqrt{1 - \alpha^2}$ . In this case $A_+^2 = A_-^2 = 1$ , so all solutions have the desired stability.

Thus to have a stable solution, we require that $\alpha^2 \leq 1$ , or $-1 \leq \alpha \leq 1$ . Plugging in the value of $\alpha$ ,

$-1 \leq \frac{2 v^2 (\Delta t)^2}{(\Delta x)^2} \sin \left(\frac{k \Delta x}{2} \right)^2 - 1 \leq 1$ $0 \leq \frac{v^2 (\Delta t)^2}{(\Delta x)^2} \sin \left(\frac{k \Delta x}{2} \right)^2 \leq 1$ $\frac{v \Delta t}{\Delta x} \left| \sin \left(\frac{k \Delta x}{2} \right) \right| \leq 1$ $v \left| \sin \left(\frac{k \Delta x}{2} \right) \right| \leq \frac{\Delta x}{\Delta t}$

In the worst case, $|\sin(k \Delta x / 2)| = 1$ , so if $v \leq \Delta x / \Delta t$ we know we’re ok for any $k$ . Interestingly, if $|\sin(k \Delta x / 2)| < 1$ for a particular mode $k$ , the physical velocity can to exceed the simulation velocity by some amount. However due to numeric error we would expect small amounts of other modes to appear at some point, which would then quickly be amplified and explode.

(e)

Do different modes decay at different rates for the stable case?

For any $k$ , the amplitude of both solutions must be exactly one. So no modes decay or grow.

(f)

Numerically solve the wave equation for the evolution from an initial condition with $u = 0$ except for one nonzero node, and verify the stability criterion.

Rearranging the relation we found in part (a),

$\begin{aligned} u_j^{n+1} &= v^2 \frac{(\Delta t)^2}{(\Delta x)^2} \left( u_{j-1}^n - 2u_j^n + u_{j+1}^n \right) +2u_j^n - u_j^{n-1} \\ &= \beta u_{j-1}^n + 2 (1 - \beta) u_j^n + \beta u_{j+1}^n - u_j^{n-1} \end{aligned}$

where we define $\beta = v^2 (\Delta t)^2 / (\Delta x)^2$ . My implementation is here.

(g)

If the equation is replaced by

$\frac{\partial^2 u}{\partial t^2} = v^2 \frac{\partial^2 u}{\partial x^2} + \gamma \frac{\partial}{\partial t} \frac{\partial^2 u}{\partial x^2}$

assume that

$u(x, t) = A e^{i(k x - \omega t)}$

and find a relationship between $k$ and $\omega$ , and simplify it for small $\gamma$ . Comment on the relationship to the preceding question.

Our ansatz is separable.

$u(x, t) = A e^{i k x} e^{-i \omega t}$

Plugging in derivatives, the differential equation becomes

$-\omega^2 = v^2 (-k^2) + \gamma (-k^2) (-i \omega) \\ -\omega^2 = -v^2 k^2 + i \gamma k^2 \omega \\ \omega^2 + i \gamma k^2 \omega -v^2 k^2 = 0 \\$

Using the quadratic formula,

$\begin{aligned} \omega &= \frac{-i \gamma k^2 \pm \sqrt{-\gamma^2 k^4 + 4 v^2 k^2}}{2} \\ &= \frac{-i \gamma k^2 \pm |k| \sqrt{4 v^2 - \gamma^2 k^2}}{2} \end{aligned}$

Manipulating to look like $\sqrt{1 + x}$ then expanding to first order in $\gamma$ ,

$\begin{aligned} \omega &= \frac{-i \gamma k^2}{2} \pm |vk| \sqrt{1 - \frac{\gamma^2 k^2}{4 v^2}} \\ &\approx \frac{-i \gamma k^2}{2} \pm |vk| \left( 1 - \frac{\gamma^2 k^2}{8 v^2} \right) \\ &\approx \frac{-i \gamma k^2}{2} \pm |vk| \end{aligned}$

This has a real part and an imaginary part. When we plug it into $u(x, t)$ the real part will cause oscillation in time, and the imaginary part will cause the amplitude to decay exponentially in time. (Note the signs; for $\gamma > 0$ it can never cause exponential growth). The damping depends on $\gamma$ , but also on $k^2$ . So the damping affects high frequencies (thus large $k$ ) much more than low ones. In the previous problem there was no damping at all.

(h)

Repeat the numerical solution of the wave equation with the same initial con- ditions, but include the damping term.

I went down a deep but ultimately very productive rabbit hole on the derivation of finite difference schemes. To cut to the chase, I now have a SymPy script that discretizes PDEs for me.

It gives two reasonable options to use for this problem. The first is an explicit scheme that is second order in space, and first order in time.

$\begin{aligned} u_{j}^{n+1} &= 2 u_{j}^{n} - u_{j}^{n-1} + v^2 \frac{(\Delta t)^2}{(\Delta x)^2} \left( u_{j-1}^{n} - 2 u_{j}^{n} + u_{j+1}^{n} \right) \\ &\phantom{=} + \gamma \frac{\Delta t}{(\Delta x)^2} \left( (u_{j-1}^{n} - 2 u_{j}^{n} + u_{j+1}^{n}) -(u_{j-1}^{n-1} - 2 u_{j}^{n-1} + u_{j+1}^{n-1}) \right) \\ &= 2 u_{j}^{n} - u_{j}^{n-1} + \left( v^2 \frac{(\Delta t)^2}{(\Delta x)^2} + \gamma \frac{\Delta t}{(\Delta x)^2} \right) \left( u_{j-1}^{n} - 2 u_{j}^{n} + u_{j+1}^{n} \right) \\ &\phantom{=} - \gamma \frac{\Delta t}{(\Delta x)^2} \left( u_{j-1}^{n-1} - 2 u_{j}^{n-1} + u_{j+1}^{n-1} \right) \end{aligned}$

The second is an implicit scheme that is second order in both space and time.

$\begin{aligned} 0 &= 2 \left(v^{2} \frac{(\Delta t)^2}{(\Delta x)^2} - 1 \right) u_{j}^{n} - v^2 \frac{(\Delta t)^2}{(\Delta x)^2} \left(u_{j-1}^{n} + u_{j+1}^{n}\right) \\ &\phantom{=} + \left(\gamma \frac{\Delta t}{(\Delta x)^2} - 1 \right) \left( u_{j}^{n+1} - u_{j}^{n-1} \right) \\ &\phantom{=} - \frac{\gamma}{2} \frac{\Delta t}{(\Delta x)^2} \left(u_{j-1}^{n+1} - u_{j-1}^{n-1} + u_{j+1}^{n+1} - u_{j+1}^{n-1}\right) \end{aligned}$

I used the explicit scheme since it’s faster to implement. I added it to my code from part (f).

8.2

Write a program to solve a 1D diffusion problem on a lattice of 500 sites, with an initial condition of zero at all the sites, except the central site which starts at the value 1.0. Take $D = \Delta x = 1$ , and use fixed boundary conditions set equal to zero.

(a)

Use the explicit finite difference scheme, and look at the behavior for $\Delta t = 1$ , 0.5, and 0.1. What step size is required by the Courant condition?

There are many finite difference schemes one might try. The standard explicit scheme for this problem uses a centered second order approximation for $x$ , and a first order backward approximation for $t$ .

$u_j^{n+1} = \left( 1 - \frac{2 D \Delta t}{(\Delta x)^2} \right) u_j^n + \frac{D \Delta t}{(\Delta x)^2} \left( u_{j-1}^n + u_{j+1}^n \right)$

The leading error terms are

$- (\Delta x)^2 \frac{D}{12} \frac{\partial^4}{\partial x^4} f(x, t) + (\Delta t)^2 \frac{1}{6} \frac{\partial^3}{\partial t^3} f(x, t) + \Delta t \frac{1}{2} \frac{\partial^2}{\partial t^2} f(x, t)$

Performing a quick stability analysis, we find that

$\begin{aligned} A &= \left( 1 - \frac{2 D \Delta t}{(\Delta x)^2} \right) + \frac{D \Delta t}{(\Delta x)^2} \left( e^{-i k \Delta x} + e^{i k \Delta x} \right) \\ &= 1 - 4 D \frac{\Delta t}{(\Delta x)^2} \sin \left( \frac{k \Delta x}{2} \right)^2 \end{aligned}$

This is between -1 and 1 when

$0 \leq D \leq \frac{1}{2} \frac{(\Delta x)^2}{\Delta t}$

(b)

Now repeat this using implicit finite differences and compare the stability.

The straightforward scheme here replaces the backward difference with a forward one (for $t$ ).

$u_j^n = \left( \frac{2 D \Delta t}{(\Delta x)^2} + 1 \right) u_j^{n+1} - \frac{D \Delta t}{(\Delta x)^2} \left( u_{j-1}^{n+1} + u_{j+1}^{n+1} \right)$

Its leading error terms are

$- (\Delta x)^2 \frac{D}{12} \frac{\partial^4}{\partial x^4} f(x, t) + (\Delta t)^2 \left( \frac{1}{6} \frac{\partial^3}{\partial t^3} f(x, t) - \frac{D}{2} \frac{\partial^4}{\partial x^2 \partial t^2} f(x, t) \right) \\ + \Delta t \left( \frac{1}{2} \frac{\partial^2}{\partial t^2} f(x, t) - D \frac{\partial^3}{\partial x^2 \partial t} f(x, t) \right)$

The important difference is that it is stable for any (positive) amount of diffusion.

$\left( \frac{2 D \Delta t}{(\Delta x)^2} + 1 \right) A - \frac{D \Delta t}{(\Delta x)^2} A \left( e^{-i k \Delta x} + e^{i k \Delta x} \right) = 1 \\ A = \frac{1}{1 + \frac{4 D \Delta t}{(\Delta x)^2} \sin \left( \frac{k \Delta x}{2} \right)^2} \\$

This is between -1 and 1 when

$0 \leq \frac{4 D \Delta t}{(\Delta x)^2} \sin \left( \frac{k \Delta x}{2} \right)^2 \\$

which reduces to

$D \geq 0$

when we assume the worst case for $\sin^2$ .

To implement this for fixed boundary conditions, we need to solve a system of linear equations that looks like this (here written for only five samples in $x$ , but the structure should be clear).

$\begin{bmatrix} 1 & 0 & 0 & 0 & 0 \\ -\alpha & 2\alpha + 1 & -\alpha & 0 & 0 \\ 0 & -\alpha & 2\alpha + 1 & -\alpha & 0 \\ 0 & 0 & -\alpha & 2\alpha + 1 & -\alpha \\ 0 & 0 & 0 & 0 & 1 \end{bmatrix} \begin{bmatrix} u_1^{n+1} \\ u_2^{n+1} \\ u_3^{n+1} \\ u_4^{n+1} \\ u_5^{n+1} \end{bmatrix} = \begin{bmatrix} u_1^n \\ u_2^n \\ u_3^n \\ u_4^n \\ u_5^n \end{bmatrix}$

Where

$\alpha = \frac{D \Delta t}{(\Delta x)^2}$

This is a specific case of the general tridiagonal system.

$\begin{bmatrix} b_1 & c_1 & 0 & 0 & 0 \\ a_2 & b_2 & c_2 & 0 & 0 \\ 0 & a_3 & b_3 & c_3 & 0 \\ 0 & 0 & a_4 & b_4 & c_4 \\ 0 & 0 & 0 & a_5 & b_5 \end{bmatrix} \begin{bmatrix} x_1 \\ x_2 \\ x_3 \\ x_4 \\ x_5 \end{bmatrix} = \begin{bmatrix} y_1 \\ y_2 \\ y_3 \\ y_4 \\ y_5 \end{bmatrix}$

We can use Gaussian elimination to systematically remove all off-diagonal elements. To start, divide the first row by $b_1$ . Then subtract $a_2$ times the first row from the second. This gets rid of $a_2$ .

$\begin{bmatrix} 1 & c_1 / b_1 & 0 & 0 & 0 \\ 0 & b_2 - a_2 \frac{c_1}{b_1} & c_2 & 0 & 0 \\ 0 & a_3 & b_3 & c_3 & 0 \\ 0 & 0 & a_4 & b_4 & c_4 \\ 0 & 0 & 0 & a_5 & b_5 \end{bmatrix} \begin{bmatrix} x_1 \\ x_2 \\ x_3 \\ x_4 \\ x_5 \end{bmatrix} = \begin{bmatrix} y_1 / b_1 \\ y_2 - a_2 \frac{y_1}{b_1} \\ y_3 \\ y_4 \\ y_5 \end{bmatrix}$

Let’s now divide the second row by $b_2 - a_2 c_1 / b_1$ , so that its leading nonzero term is one.

$\begin{bmatrix} 1 & c_1 / b_1 & 0 & 0 & 0 \\ 0 & 1 & \frac{c_2}{b_2 - a_2 c_1 / b_1} & 0 & 0 \\ 0 & a_3 & b_3 & c_3 & 0 \\ 0 & 0 & a_4 & b_4 & c_4 \\ 0 & 0 & 0 & a_5 & b_5 \end{bmatrix} \begin{bmatrix} x_1 \\ x_2 \\ x_3 \\ x_4 \\ x_5 \end{bmatrix} = \begin{bmatrix} y_1 / b_1 \\ \frac{y_2 - a_2 y_1 / b_1}{b_2 - a_2 c_1 / b_1} \\ y_3 \\ y_4 \\ y_5 \end{bmatrix}$

We can repeat this process until all entries below the diagonal are eliminated, and all diagonal entries are one. The result is

$\begin{bmatrix} 1 & c_1' & 0 & 0 & 0 \\ 0 & 1 & c_2' & 0 & 0 \\ 0 & 0 & 1 & c_3' & 0 \\ 0 & 0 & 0 & 1 & c_4' \\ 0 & 0 & 0 & 0 & 1 \end{bmatrix} \begin{bmatrix} x_1 \\ x_2 \\ x_3 \\ x_4 \\ x_5 \end{bmatrix} = \begin{bmatrix} y_1' \\ y_2' \\ y_3' \\ y_4' \\ y_5' \end{bmatrix}$

where

$\begin{aligned} c_i' &= \begin{cases} \Large{\frac{c_1}{b_1}} &\text{for } i = 1 \\ \Large{\frac{c_i}{b_i - a_i c_{i-1}'}} &\text{for } 1 < i < m \end{cases} \\ y_i' &= \begin{cases} \Large{\frac{y_1}{b_1}} &\text{for } i = 1 \\ \Large{\frac{y_i - a_i y_{i - 1}'}{b_i - a_i c_{i-1}'}} &\text{for } 1 < i \leq m \end{cases} \\ \end{aligned}$

Now we remove the entries above the diagonal, starting with $c_4'$ . Multiply the last row by $c_4'$ and subtract from the fourth. Note that this doesn’t change any diagonal entries. After doing this $m - 1$ times we have

$\begin{bmatrix} 1 & 0 & 0 & 0 & 0 \\ 0 & 1 & 0 & 0 & 0 \\ 0 & 0 & 1 & 0 & 0 \\ 0 & 0 & 0 & 1 & 0 \\ 0 & 0 & 0 & 0 & 1 \end{bmatrix} \begin{bmatrix} x_1 \\ x_2 \\ x_3 \\ x_4 \\ x_5 \end{bmatrix} = \begin{bmatrix} y_1'' \\ y_2'' \\ y_3'' \\ y_4'' \\ y_5'' \end{bmatrix}$

where

$y_i'' = \begin{cases} y_m' &\text{for } i = m \\ y_i' - c_i' y_{i+1}'' &\text{for } 1 \leq i < m \end{cases} \\$

These $y_i''$ are the solutions we wanted.

For our particular system,

$\begin{aligned} a_i &= \begin{cases} -\alpha &\text{for } 1 < i < m \\ 0 &\text{for } i = m \end{cases} \\ b_i &= \begin{cases} 1 &\text{for } i = 1 \\ 2 \alpha + 1 &\text{for } 1 < i < m \\ 1 &\text{for } i = m \end{cases} \\ c_i &= \begin{cases} 0 &\text{for } i = 1 \\ -\alpha &\text{for } 1 < i < m \end{cases} \\ x_i &= u_i^{n + 1} \\ y_i &= u_i^n \end{aligned}$

We can compute the $c_i'$ in advance and reuse them. Then $y_i'$ and $y_i''$ reduce to

$\begin{aligned} y_i' &= \begin{cases} y_1 &\text{for } i = 1 \\ \Large{\frac{y_i + \alpha y_{i-1}'}{1 + \alpha (2 + c_i')}} &\text{for } 1 < i < m \\ y_n &\text{for } i = m \end{cases} \\ \end{aligned} \\ \begin{aligned} y_i'' &= \begin{cases} y_m &\text{for } i = m \\ y_i' - c_i' y_{i+1}'' &\text{for } m > i > 1 \\ y_1 &\text{for } i = 1 \end{cases} \\ \end{aligned}$

Interestingly, by using a slightly different Gaussian elimination scheme, we can get rid of the auxiliary array of $c_i'$ values. I derive this in my notes.