Finite Differences: PDEs

Problem 8.1: Wave Equation

Consider the wave equation

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

We can approximate this with a finite difference equation:

$$ \frac{(u^{n+1}_j - u^{n}_j) - (u^{n}_j - u^{n-1}_j)}{(\Delta t)^2} = v^2 \frac{(u^{n}_{j+1} - u^{n}_{j}) - (u^{n}_{j} - u^{n}_{j-1})}{(\Delta x)^2} $$

Solving for $u^{n+1}_j$, we get:

$$u^{n+1} = 2u^n_j - u^{n-1}_j + \Big(\frac{v \Delta t}{\Delta x}\Big)^2 [u^n_{j+1} - 2 u^n_j + u^n_{j-1}]$$

Consider ansatz $u^{n}_j = A(k)^n e^{ijk \Delta x}$. Plugging in the previous equation, we get:

$$A^{n+1} e^{ijk \Delta x} = 2 A^n e^{ijk \Delta x} - A^{n-1} e^{ijk \Delta x} - \Big(\frac{v \Delta t}{\Delta x}\Big)^2 [ A^n e^{ik (j+1) \Delta x} - 2 A^n e^{ijk \Delta x} + A^n e^{ik (j-1) \Delta x}] $$
$$A e^{ijk \Delta x} = 2 e^{ijk \Delta x} - \frac{1}{A} e^{ijk \Delta x} - \Big(\frac{v \Delta t}{\Delta x}\Big)^2 [ e^{ik (j+1) \Delta x} - 2 e^{ijk \Delta x} + e^{ik (j-1) \Delta x}] $$
$$A = 2 - \frac{1}{A} - \Big(\frac{v \Delta t}{\Delta x}\Big)^2 [ e^{ik \Delta x} - 2+ e^{-ik \Delta x}] $$
$$A = 2 - \frac{1}{A} - \Big(\frac{v \Delta t}{\Delta x}\Big)^2 [- 2 + 2 \cos(k \Delta x)] $$
$$A = 2 - \frac{1}{A} - \omega$$

Where we substituted $\omega = \Big(\frac{v \Delta t}{\Delta x}\Big)^2 [- 2 + 2 \cos(k \Delta x)] $

$$A^2 = 2 A- 1- A \omega = A (2 - \omega) - 1$$
$$A = \frac{1}{2} (2 \pm \sqrt{-4 \omega+ \omega^2} - \omega) $$

There are two different integer solutions, with $A = -1, ω = 4 $ and $A = 1, ω = 0$

Now to find conditions on velocity time step and space step for stability, we first note that $$A_0 A_1= \frac{1}{4} (2 + \sqrt{-4 + \omega} \sqrt{\omega} - \omega) (2 - \sqrt{-4 + \omega} \sqrt{\omega} - \omega) = 1$$

Which means that $A_0$ and $A_1$ must have magnitude $1$, since otherwise they would either decay or explode.

We note that different modes don't decay at different time rates, since the only dependence on $j$, the time factor, is though $\cos(k \Delta x)$, which doesn't decay.

$$u^{n+1}_j = 2u^n_j - u^{n-1}_j + \Big(\frac{v \Delta t}{\Delta x}\Big)^2 [u^n_{j+1} - 2 u^n_j + u^n_{j-1}]$$
In [ ]:
# Solve numerically the wave equation

import matplotlib.pyplot as plt
%matplotlib notebook

import numpy as np

fig = plt.figure()
ax = fig.add_subplot(111)
fig.show()
fig.canvas.draw()

v = 1
dt = 0.1
dx = 0.2
N = 500
H = 10
T = 10000

x = np.arange(0, N)
past = np.zeros(N)
present = np.zeros(N)
future = np.zeros(N)

present[N//2] = H - 2

α = (v * dt / dx) ** 2  
    
for t in range(T):
    future = 2 * present - past
    for j in range(1, N-1):
        future[j] += α * (present[j + 1] - 2 * present[j] + present[j-1])
    past, present = present, future
    
    if t % 20 == 0:
        ax.clear()
        ax.plot(x, present)
        fig.canvas.draw()

Now we replace the equation 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 now $u(x,t) = A e^{i(kx - \omega t)}$

In [2]:
from scipy.ndimage.filters import convolve

%matplotlib notebook
fig = plt.figure()
ax = fig.add_subplot(111)

fig.show()
fig.canvas.draw()

N = 100
T = 10000

U = np.zeros([N, N])
U[-1, :] = -1
U[:, -1] = 1

V = np.copy(U)
plt.imshow(V)
plt.show()

α = 1.5

for t in range(T):
    for k in range(1, N-1):
        for j in range(1, N-1):
            V[j][k] = (1-α) * U[j][k]
            
            V[j][k] += α/4 * U[j][k+1]
            V[j][k] += α/4 * U[j+1][k]
            
            V[j][k] += α/4 * V[j][k-1]
            V[j][k] += α/4 * V[j-1][k]

    
    U = np.copy(V)
    
    if t % 2 == 0:
        ax.clear()
        ax.imshow(V, cmap='hot')
        fig.canvas.draw()

---------------------------------------------------------------------------
KeyboardInterrupt                         Traceback (most recent call last)
<ipython-input-2-ccded344605e> in <module>
     26             V[j][k] = (1-α) * U[j][k]
     27 
---> 28             V[j][k] += α/4 * U[j][k+1]
     29             V[j][k] += α/4 * U[j+1][k]
     30 

KeyboardInterrupt:
In [ ]: