MATMEK-4270 presentations – Finite difference methods for the wave equation

The wave equation is a partial differential equation (PDE)

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

We will consider the time and space domains: $t \in [0, T]$ , $x \in [0, L]$ .

The wave equation is an initial-boundary value problem!
Two initial conditions required since two derivatives in time
Two boundary conditions required since two derivatives in space
The solutions are waves that can be written as $u (x + c t)$ and $u (x - c t)$

Wave solution with different boundary conditions

Boundary conditions

Dirichlet (Fixed end)

$u (0, t) = u (L, t) = 0$

The wave will be reflected, but $u$ will change sign. A nonzero Dirichlet condition is also possible, but will not be considered here.

Neumann (Loose end)

$\frac{\partial u}{\partial x} (0, t) = \frac{\partial u}{\partial x} (L, t) = 0$

The wave will be reflected without change in sign. A nonzero Neumann condition is also possible, but will not be considered here.

Boundary conditions continued

Open boundary (No end)

$\frac{\partial u (0, t)}{\partial t} - c \frac{\partial u (0, t)}{\partial x} = 0$

$\frac{\partial u (L, t)}{\partial t} + c \frac{\partial u (L, t)}{\partial x} = 0$

The wave will simply pass undisturbed and unreflected through an open boundary.

Periodic boundary (No end)

$u (0, t) = u (L, t)$

The solution repeats itself indefinitely.

Discretization

The simplest possible discretization is uniform in time and space

$t_{n} = n Δ t, n = 0, 1, \dots, N_{t}$

$x_{j} = j Δ x, j = 0, 1, \dots, N$

A mesh function in space and time is defined as

$u_{j}^{n} = u (x_{j}, t_{n})$

The mesh function has one value at each node in the mesh. For simplicity in later algorithms we will use the vectors

$u^{n} = (u_{0}^{n}, u_{1}^{n}, \dots, u_{N}^{n})^{T},$

which is the solution vector at time $t_{n}$ .

A second order accurate discretization of the wave equation is

$\frac{u_{j}^{n + 1} - 2 u_{j}^{n} + u_{j}^{n - 1}}{Δ t^{2}} = c^{2} \frac{u_{j + 1}^{n} - 2 u_{j}^{n} + u_{j - 1}^{n}}{Δ x^{2}}$

The finite difference stencil makes use of 5 neighboring points

$\frac{u_{j}^{n + 1} - 2 u_{j}^{n} + u_{j}^{n - 1}}{Δ t^{2}} = c^{2} \frac{u_{j + 1}^{n} - 2 u_{j}^{n} + u_{j - 1}^{n}}{Δ x^{2}}$

Can only be used for internal points

The finite difference stencil is not used at the spatial boundary

$\frac{u_{N}^{n + 1} - 2 u_{N}^{n} + u_{N}^{n - 1}}{Δ t^{2}} = c^{2} \frac{u_{N + 1}^{n} - 2 u_{N}^{n} + u_{N - 1}^{n}}{Δ x^{2}}$

Used at the boundary the regular stencil will contain a ghost node
But at the boundary we use boundary conditions and do not solve the PDE!

We use a marching method in time

Initialize $u^{0}$ and $u^{1}$
for n in range( $1, N_{t} - 1$ ):
- for j in range( $1, N - 1$ ):
  - $u_{j}^{n + 1} = 2 u_{j}^{n} - u_{j}^{n - 1} + {(\frac{c Δ t}{Δ x})}^{2} (u_{j + 1}^{n} - 2 u_{j}^{n} + u_{j - 1}^{n})$

and apply the chosen boundary conditions.

All the indices makes it a bit messy. Lets make use of a differentiation matrix for the spatial dimension! And the Courant (or CFL) number

$\overset{―}{c} = \frac{c Δ t}{Δ x}$

Use differentiation matrix to simplify the notation

We define the second differentiation matrix without the scaling $1 / (Δ x)^{2}$ such that

$D^{(2)} = [\begin{matrix} - 2 & 1 & 0 & 0 & 0 & 0 & 0 & 0 \\ 1 & - 2 & 1 & 0 & 0 & 0 & 0 & \dots \\ 0 & 1 & - 2 & 1 & 0 & 0 & 0 & \dots \\ ⋮ & ⋱ & \dots \\ ⋮ & 0 & 0 & 0 & 1 & - 2 & 1 & 0 \\ ⋮ & 0 & 0 & 0 & 0 & 1 & - 2 & 1 \\ 0 & 0 & 0 & 0 & 0 & 0 & 1 & - 2 \end{matrix}]$

and thus row $0 < j < N$ of $D^{(2)} u^{n}$ becomes

$(D^{(2)} u^{n})_{j} = u_{j + 1}^{n} - 2 u_{j}^{n} + u_{j - 1}^{n}$

The vectorized marching method becomes

Initialize $u^{0}$ and $u^{1}$
for n in range( $1, N_{t} - 1$ ):
- $u^{n + 1} = 2 u^{n} - u^{n - 1} + {\underset{―}{c}}^{2} D^{(2)} u^{n}$
- Apply boundary conditions to $u_{0}^{n + 1}$ and $u_{N}^{n + 1}$

The boundary step can often, but not always, be incorporated into the matrix $D^{(2)}$
Very easy to vectorize using the matrix vector product!

PDE solvers (of time-dependent problems) should use memory carefully

Note

At any time we only need to store three vectors: $u^{n + 1}, u^{n}$ and $u^{n - 1}$ .
- Memory requirement = $3 (N + 1)$ floating point numbers
Storing all time steps requires $(N_{t} + 1) \times (N + 1)$ floating point numbers
Not a huge problem for our case, but for 2 or 3 spatial dimensions it is very important!

Implementation - A low-memory marching method needs to update solution vectors

Allocate three vectors $u^{n m 1}, u^{n}, u^{n p 1}$ , representing $u^{n - 1}, u^{n}, u^{n + 1}$ .
Initialize $u^{0}$ and $u^{1}$ by setting $u^{n m 1} = u^{0}, u^{n} = u^{1}$
for n in range( $1, N_{t} - 1$ ):
- $u^{n p 1} = 2 u^{n} - u^{n m 1} + {\underset{―}{c}}^{2} D^{(2)} u^{n}$
- Apply boundary conditions to $u_{0}^{n p 1}$ and $u_{N}^{n p 1}$
- Update to next iteration:
  - $u^{n m 1} \leftarrow u^{n}$
  - $u^{n} \leftarrow u^{n p 1}$

In Python

Set up solver

import numpy as np 
from scipy import sparse
import sympy as sp 
x, t = sp.symbols('x,t')
N = 100
Nt = 500 
L = 2
c = 1 # wavespeed
dx = L / N
CFL = 1.0
dt = CFL*dx/c
xj = np.linspace(0, L, N+1)
unm1, un, unp1 = np.zeros((3, N+1))
D2 = sparse.diags([1, -2, 1], [-1, 0, 1], (N+1, N+1))
u0 = sp.exp(-200*(x-L/2+t)**2)

Solve by marching method

unm1[:] = sp.lambdify(x, u0.subs(t, 0))(xj)
un[:] = sp.lambdify(x, u0.subs(t, dt))(xj) 
for n in range(Nt):
  unp1[:] = 2*un - unm1 + CFL**2 * D2 @ un 
  unp1[0] = 0
  unp1[-1] = 0
  unm1[:] = un 
  un[:] = unp1

Store results at intermediate intervals for plotting

unm1[:] = sp.lambdify(x, u0.subs(t, 0))(xj)
un[:] = sp.lambdify(x, u0.subs(t, dt))(xj) 
plotdata = {0: unm1.copy()}
for n in range(Nt):
  unp1[:] = 2*un - unm1 + CFL**2 * D2 @ un 
  unp1[0] = 0
  unp1[-1] = 0
  unm1[:] = un 
  un[:] = unp1
  if n % 10 == 0:
    plotdata[n] = unp1.copy()

For example every tenth time step. Normally you do not need every time step to get a good animation.

Create animation after the simulation is finished

Animated Image

def animation(data):
  from matplotlib import animation 
  fig, ax = plt.subplots()
  v = np.array(list(data.values()))
  t = np.array(list(data.keys()))
  save_step = t[1]-t[0]
  line, = ax.plot(xj, data[0])
  ax.set_ylim(v.min(), v.max())
  def update(frame):
    line.set_ydata(data[frame*save_step])
    return (line,)
  ani = animation.FuncAnimation(fig=fig, func=update, frames=len(data), blit=True)
  ani.save('wavemovie.apng', writer='pillow', fps=5) # This animated png opens in a browser

How to implement the initial conditions?

To initialize a mesh function $u^{0}$ , we write $u^{0} = I (x)$

which represents

$u_{j}^{0} = I (x_{j}), \forall j = 0, 1, \dots, N$

u0 = sp.exp(-200*(x-L/2+t)**2) 
unm1[:] = sp.lambdify(x, u0.subs(t, 0))(xj)

How about the second condition $\frac{\partial u}{\partial t} (x, 0) = 0$ ?

Just like for the vibration equation there are several options. If you have an analytical solution $I (x, t)$ that is time-dependent you can specify one wave as:

$u^{1} = I (x, Δ t)$

un[:] = sp.lambdify(x, u0.subs(t, dt))(xj)

If you do not have $I (x, t)$ , then what?

How to fix $\frac{\partial u}{\partial t} (x, 0) = 0$ , option 1

Use a forward difference

$\frac{\partial u}{\partial t} (x, 0) \approx \frac{u^{1} - u^{0}}{Δ t} = 0, such that u^{1} = u^{0}$

Only first order accurate, but still a possibility.

Use a second order forward difference

$\frac{\partial u}{\partial t} (x, 0) \approx \frac{- u^{2} + 4 u^{1} - 3 u^{0}}{2 Δ t} = 0, such that u^{1} = \frac{3 u^{0} + u^{2}}{4}$

Second order accurate, but implicit. And you get two waves!

How to implement $\frac{\partial u}{\partial t} (x, 0) = 0$ , option 2

Use a second order central differece

$\frac{\partial u}{\partial t} (x, 0) = \frac{u^{1} - u^{- 1}}{2 Δ t} = 0, such that u^{1} = u^{- 1}$

and the PDE at $n = 0$

$u^{1} = 2 u^{0} - u^{- 1} + {\underset{―}{c}}^{2} D^{(2)} u^{0}$

Insert for $u^{- 1} = u^{1}$ to obtain

$u^{1} = u^{0} + \frac{{\underset{―}{c}}^{2}}{2} D^{(2)} u^{0}$

Second order accurate and explicit

How to fix boundary conditions?

We will consider 4 different types of boundary conditions

Dirichlet	$u (0, t)$ and $u (L, t)$
Neumann	$\frac{\partial u}{\partial x} (0, t)$ and $\frac{\partial u}{\partial x} (L, t)$
Open	$\frac{\partial u}{\partial t} (0, t) - c \frac{\partial u}{\partial x} (0, t) = 0$ and $\frac{\partial u}{\partial t} (L, t) + c \frac{\partial u}{\partial x} (L, t) = 0$
Periodic	$u (L, t) = u (0, t)$

Note

Accounting for boundary conditions very often takes more than 50 % of the lines of code in a PDE solver!

Dirichlet boundary conditions

We need to fix $u (0, t) = I (0)$ and $u (L, t) = I (L)$ and start by fixing this at $t = 0$

$u_{0}^{0} = I (0) and u_{N}^{0} = I (L)$

Next, we compute

$u^{1} = u^{0} + \frac{{\underset{―}{c}}^{2}}{2} D^{(2)} u^{0}$

Here, if the first and last rows of $D^{(2)}$ are set to zero, then $u_{0}^{1} = u_{0}^{0}$ and $u_{N}^{1} = u_{N}^{0}$ .

Next, for $n = 1, 2, \dots, N_{t} - 1$

$u^{n + 1} = 2 u^{n} - u^{n - 1} + {\underset{―}{c}}^{2} D^{(2)} u^{n}$

Again, if the first and last rows of $D^{(2)}$ are zero, then $u_{0}^{n + 1} = u_{0}^{0}$ and $u_{N}^{n + 1} = u_{N}^{0}$ for all $n$ . The boundary values remain as initially set at $t = 0$ .

Dirichlet boundary conditions summary

Set $u^{0} = I (x)$ and define a modified differentiation matrix

${\tilde{D}}^{(2)} = [\begin{matrix} 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 \\ 1 & - 2 & 1 & 0 & 0 & 0 & 0 & \dots \\ 0 & 1 & - 2 & 1 & 0 & 0 & 0 & \dots \\ ⋮ & ⋱ & \dots \\ ⋮ & 0 & 0 & 0 & 1 & - 2 & 1 & 0 \\ ⋮ & 0 & 0 & 0 & 0 & 1 & - 2 & 1 \\ 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 \end{matrix}]$

Now boundary conditions will be ok at all time steps simply by:

Initialize $u^{0}$ and compute $u^{1} = u^{0} + \frac{{\underset{―}{c}}^{2}}{2} {\tilde{D}}^{(2)} u^{0}$ . Set $u^{n m 1} = u^{0}, u^{n} = u^{1}$
for n in range( $1, N_{t} - 1$ ):
- $u^{n p 1} = 2 u^{n} - u^{n m 1} + {\underset{―}{c}}^{2} {\tilde{D}}^{(2)} u^{n}$
- Update to next iteration: $u^{n m 1} = u^{n}; u^{n} = u^{n p 1}$

Dirichlet boundary conditions summary

Note

It is also possible to do nothing with $D^{(2)}$ and simply fix the boundary conditions after updating all the internal points

Initialize $u^{n m 1} = u^{0}$ and compute $u^{n} = u^{1} = u^{0} + \frac{{\underset{―}{c}}^{2}}{2} D^{(2)} u^{0}$ .
Set $u_{0}^{n m 1} = u_{0}^{n} = 0$ and $u_{N}^{n m 1} = u_{N}^{n} = 0$ .
for n in range( $1, N_{t} - 1$ ):
- $u^{n p 1} = 2 u^{n} - u^{n m 1} + {\underset{―}{c}}^{2} D^{(2)} u^{n}$
- Set $u_{0}^{n p 1} = 0$ and $u_{N}^{n p 1} = 0$
- Update to next iteration: $u^{n m 1} \leftarrow u^{n}; u^{n} \leftarrow u^{n p 1}$

Note

Regular, unmodified $D^{(2)}$ , where the first and last rows are completely irrelevant.

Neumann boundary conditions

We need to fix $\frac{\partial u}{\partial x} (0, t) = 0$ and $\frac{\partial u}{\partial x} (L, t) = 0$ . We already have $u^{0} = I (x)$ .

A second order central scheme at $x = 0$ is using ghost cell at $j = - 1$

$\frac{\partial u}{\partial x} (0, t_{n}) = \frac{u_{1}^{n} - u_{- 1}^{n}}{2 Δ x} = 0 \to u_{- 1}^{n} = u_{1}^{n}$

Use the ghost cell and the PDE to fix the Neumann condition

The PDE at the left hand side $j = 0$ using ghost cell:

$u_{0}^{n + 1} = 2 u_{0}^{n} - u_{0}^{n - 1} + {\underset{―}{c}}^{2} (u_{1}^{n} - 2 u_{0}^{n} + u_{- 1}^{n})$

Insert for $u_{- 1}^{n} = u_{1}^{n}$ and obtain

$u_{0}^{n + 1} = 2 u_{0}^{n} - u_{0}^{n - 1} + {\underset{―}{c}}^{2} (2 u_{1}^{n} - 2 u_{0}^{n})$

Note

Second order accurate and explicit. Can be implemented by modifying $D^{(2)}$ !

Neumann at $x = L$ is the same

$\frac{\partial u}{\partial x} (L, t_{n}) = \frac{u_{N + 1}^{n} - u_{N - 1}^{n}}{2 Δ x} = 0 \to u_{N + 1}^{n} = u_{N - 1}^{n}$

The PDE at the right hand side $j = N$ using ghost cell:

$u_{N}^{n + 1} = 2 u_{N}^{n} - u_{N}^{n - 1} + {\underset{―}{c}}^{2} (u_{N + 1}^{n} - 2 u_{N}^{n} + u_{N - 1}^{n})$

Insert for $u_{N + 1}^{n} = u_{N - 1}^{n}$ and obtain

$u_{N}^{n + 1} = 2 u_{N}^{n} - u_{N}^{n - 1} + {\underset{―}{c}}^{2} (2 u_{N - 1}^{n} - 2 u_{N}^{n})$

And for $n = 1$ we similarly get

$u_{0}^{1} = u_{0}^{0} + \frac{{\underset{―}{c}}^{2}}{2} (2 u_{1}^{n} - 2 u_{0}^{n}) and u_{N}^{1} = u_{N}^{0} + \frac{{\underset{―}{c}}^{2}}{2} (2 u_{N - 1}^{n} - 2 u_{N}^{n})$

Neumann summary

Set $u^{0} = I (x)$ and define a modified differentiation matrix

${\tilde{D}}^{(2)} = [\begin{matrix} - 2 & 2 & 0 & 0 & 0 & 0 & 0 & 0 \\ 1 & - 2 & 1 & 0 & 0 & 0 & 0 & \dots \\ 0 & 1 & - 2 & 1 & 0 & 0 & 0 & \dots \\ ⋮ & ⋱ & \dots \\ ⋮ & 0 & 0 & 0 & 1 & - 2 & 1 & 0 \\ ⋮ & 0 & 0 & 0 & 0 & 1 & - 2 & 1 \\ 0 & 0 & 0 & 0 & 0 & 0 & 2 & - 2 \end{matrix}]$

Now boundary conditions will be ok at all time steps simply by:

Initialize $u^{0}$ and compute $u^{1} = u^{0} + \frac{{\underset{―}{c}}^{2}}{2} {\tilde{D}}^{(2)} u^{0}$ . Set $u^{n m 1} = u^{0}, u^{n} = u^{1}$
for n in range( $1, N_{t} - 1$ ):
- $u^{n p 1} = 2 u^{n} - u^{n m 1} + {\underset{―}{c}}^{2} {\tilde{D}}^{(2)} u^{n}$
- Update to next iteration: $u^{n m 1} \leftarrow u^{n}; u^{n} \leftarrow u^{n p 1}$

Open boundary

The wave simply disappears through the boundary

$\frac{\partial u}{\partial t} (0, t) - c \frac{\partial u}{\partial x} (0, t) = 0 and \frac{\partial u}{\partial t} (L, t) + c \frac{\partial u}{\partial x} (L, t) = 0$

As for Neumann there are several ways to implement these boundary conditions. The simplest option is to solve the first order accurate

$\frac{u_{0}^{n + 1} - u_{0}^{n}}{Δ t} - c \frac{u_{1}^{n} - u_{0}^{n}}{Δ x} = 0$

such that

$u_{0}^{n + 1} = u_{0}^{n} + \frac{c Δ t}{Δ x} (u_{1}^{n} - u_{0}^{n})$

Second order option

$\frac{u_{0}^{n + 1} - u_{0}^{n - 1}}{2 Δ t} - c \frac{- u_{2}^{n} + 4 u_{1}^{n} - 3 u_{0}^{n}}{2 Δ x} = 0$

Solve for the boundary node $u_{0}^{n + 1}$

$u_{0}^{n + 1} = u_{0}^{n - 1} + \frac{c Δ t}{Δ x} (- u_{2}^{n} + 4 u_{1}^{n} - 3 u_{0}^{n})$

Nice option, but difficult to incorporate in the $D^{(2)}$ matrix, since there is no way to modify the first and last rows of $D^{(2)}$ such that

$u_{0}^{n + 1} = 2 u_{0}^{n} - u_{0}^{n - 1} + {\underset{―}{c}}^{2} (D^{(2)} u^{n})_{0}$

Second second order option

Use central, second order scheme

$\frac{u_{0}^{n + 1} - u_{0}^{n - 1}}{2 Δ t} - c \frac{u_{1}^{n} - u_{- 1}^{n}}{2 Δ x} = 0$

and isolate the ghost node $u_{- 1}^{n}$ :

$u_{- 1}^{n} = u_{1}^{n} - \frac{1}{\underset{―}{c}} (u_{0}^{n + 1} - u_{0}^{n - 1})$

Use regular PDE at the boundary that includes the ghost node:

$u_{0}^{n + 1} = 2 u_{0}^{n} - u_{0}^{n - 1} + {\underset{―}{c}}^{2} (u_{1}^{n} - 2 u_{0}^{n} + u_{- 1}^{n})$

This gives an equation for $u_{0}^{n + 1}$ that fixes the open boundary condition:

$u_{0}^{n + 1} = 2 (1 - \underset{―}{c}) u_{0}^{n} - \frac{1 - \underset{―}{c}}{1 + \underset{―}{c}} u_{0}^{n - 1} + \frac{2 {\underset{―}{c}}^{2}}{1 + \underset{―}{c}} u_{1}^{n_{1}}$

Open boundary conditions

Left boundary:

$u_{0}^{n + 1} = 2 (1 - \underset{―}{c}) u_{0}^{n} - \frac{1 - \underset{―}{c}}{1 + \underset{―}{c}} u_{0}^{n - 1} + \frac{2 {\underset{―}{c}}^{2}}{1 + \underset{―}{c}} u_{1}^{n_{1}}$

Right boundary:

$u_{N}^{n + 1} = 2 (1 - \underset{―}{c}) u_{N}^{n} - \frac{1 - \underset{―}{c}}{1 + \underset{―}{c}} u_{N}^{n - 1} + \frac{2 {\underset{―}{c}}^{2}}{1 + \underset{―}{c}} u_{N - 1}^{n_{1}}$

Both explicit and second order. But not possible to implement into the matrix such that

$u^{n + 1} = 2 u^{n} - u^{n - 1} + {\underset{―}{c}}^{2} D^{(2)} u^{n}$

Implementation open boundaries

Initialize $u^{0}$ and compute $u^{1} = u^{0} + \frac{{\underset{―}{c}}^{2}}{2} D^{(2)} u^{0}$ . Set $u^{n m 1} = u^{0}, u^{n} = u^{1}$
for n in range( $1, N_{t} - 1$ ):
- $u^{n p 1} = 2 u^{n} - u^{n m 1} + {\underset{―}{c}}^{2} D^{(2)} u^{n}$
- $u_{0}^{n p 1} = 2 (1 - \underset{―}{c}) u_{0}^{n} - \frac{1 - \underset{―}{c}}{1 + \underset{―}{c}} u_{0}^{n m 1} + \frac{2 {\underset{―}{c}}^{2}}{1 + \underset{―}{c}} u_{1}^{n m 1}$
- $u_{N}^{n p 1} = 2 (1 - \underset{―}{c}) u_{N}^{n} - \frac{1 - \underset{―}{c}}{1 + \underset{―}{c}} u_{N}^{n m 1} + \frac{2 {\underset{―}{c}}^{2}}{1 + \underset{―}{c}} u_{N - 1}^{n m 1}$
- Update to next iteration: $u^{n m 1} \leftarrow u^{n}; u^{n} \leftarrow u^{n p 1}$

Note

There is no need to use a modified $D^{(2)}$ . The two updates of $u_{0}^{n p 1}$ and $u_{N}^{n p 1}$ will overwrite anything computed in the first step.

Periodic boundary conditions

A periodic solution is a solution that is repeating itself indefinitely. For example $u (x) = \sin (2 π x)$ :

We solve the problem for example for $x \in [0, 1]$ , but the actual solution will be like above, with no boundaries.

Note

A periodic domain is also referred to as a domain with no boundaries.

A periodic mesh in time

$u (t_{N}) = u (0) or u^{N} = u^{0}$

Note

There are only $N$ unknowns $u^{0}, u^{1}, \dots, u^{N - 1}$ for a mesh with $N + 1$ nodes.

Consider the discretization of $u^{″}$

At the left hand side of the domain, the point to the left of $u^{0}$ is $u^{N - 1}$

$u^{″} (0) \approx \frac{u^{1} - 2 u^{0} + u^{- 1}}{h^{2}} = \frac{u^{1} - 2 u^{0} + u^{N - 1}}{h^{2}}$

At the right hand side of the domain the point to the right of $u^{N - 1}$ is $u^{N} = u^{0}$

$u^{″} (t_{N - 1}) \approx \frac{u^{N} - 2 u^{N - 1} + u^{N - 2}}{h^{2}} = \frac{u^{0} - 2 u^{N - 1} + u^{N - 2}}{h^{2}}$

Periodic boundary conditions can be implemented in the matrix $D^{(2)} \in R^{N + 1 \times N + 1}$

${\tilde{D}}^{(2)} = [\begin{matrix} - 2 & 1 & 0 & 0 & 0 & 0 & 1 & 0 \\ 1 & - 2 & 1 & 0 & 0 & 0 & 0 & \dots \\ 0 & 1 & - 2 & 1 & 0 & 0 & 0 & \dots \\ ⋮ & ⋱ & \dots \\ ⋮ & 0 & 0 & 0 & 1 & - 2 & 1 & 0 \\ ⋮ & 0 & 0 & 0 & 0 & 1 & - 2 & 1 \\ 0 & 0 & 0 & 0 & 0 & 0 & 1 & - 2 \end{matrix}]$

Note that the matrix expects $u = (u^{0}, u^{1}, \dots, u^{N - 1}, u^{N})$ , even though $u^{0} = u^{N}$ . The last row in ${\tilde{D}}^{(2)}$ is thus irrelevant, because we wil set $u^{0} = u^{N}$ manually.

Implementation periodic boundaries

Initialize $u^{0}$ and compute $u^{1} = u^{0} + \frac{{\underset{―}{c}}^{2}}{2} {\tilde{D}}^{(2)} u^{0}$ . Set $u^{n m 1} = u^{0}, u^{n} = u^{1}$
for n in range(1, $N_{t} - 1$ ):
- $u^{n p 1} = 2 u^{n} - u^{n m 1} + {\underset{―}{c}}^{2} {\tilde{D}}^{(2)} u^{n}$
- $u_{N}^{n p 1} = u_{0}^{n p 1}$
- Update to next iteration: $u^{n m 1} = u^{n}; u^{n} = u^{n p 1}$

Periodic wave

Properties of the wave equation

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

If the initial condition is $u (x, 0) = I (x)$ and $\frac{\partial u}{\partial t} (x, 0) = 0$ , then the solution at $t > 0$ is

$u (x, t) = \frac{1}{2} (I (x - c t) + I (x + c t))$

These are two waves - one traveling to the left and the other traveling to the right

If the initial condition $I (x) = e^{i k x}$ , then

$u (x, t) = \frac{1}{2} (e^{i k (x - c t)} + e^{i k (x + c t)})$

is a solution

Representation of waves as complex exponentials

If the initial condition is a sum of waves (superposition, each wave is a solution of the wave equation)

$I (x) = \sum_{k = 0}^{K} a_{k} e^{i k x} = \sum_{k = 0}^{K} a_{k} (\cos k x + i \sin k x)$

for some $K$ , then the solution is

$u (x, t) = \frac{1}{2} \sum_{k = 0}^{K} a_{k} (e^{i k (x - c t)} + e^{i k (x + c t)})$

We will analyze one component $e^{i k (x + c t)} = e^{i k x + ω t}$ , where $ω = k c$ is the frequency in time. This is very similar to the investigation we did for the numerical frequency for the vibration equation.

Assume that the numerical solution is a complex wave

$u (x_{j}, t_{n}) = u_{j}^{n} = e^{i k (x_{j} + \tilde{ω} t_{n})}$

How accurate is $\tilde{ω}$ compared to the exact $ω = k c$ ?
What can be concluded about stability?

Note that the solution is a recurrence relation

$u_{j}^{n} = e^{i k x_{j}} e^{i \tilde{ω} n Δ t} = (e^{i \tilde{ω} Δ t})^{n} e^{i k x_{j}}$

with an amplification factor $A = e^{i \tilde{ω} Δ t}$ such that

$u_{j}^{n} = A^{n} e^{i k x_{j}}$

Numerical dispersion relation

We can find $\tilde{ω}$ by inserting for $e^{i k (x_{j} + \tilde{ω} t_{n})}$ in the discretized wave equation

$\frac{u_{j}^{n + 1} - 2 u_{j}^{n} + u_{j}^{n - 1}}{Δ t^{2}} = c^{2} \frac{u_{j + 1}^{n} - 2 u_{j}^{n} + u_{j - 1}^{n}}{Δ x^{2}}$

This is a lot of work, just like it was for the vibration equation. In the end we should get

$\tilde{ω} = \frac{2}{Δ t} \sin^{- 1} (C \sin (\frac{k Δ x}{2}))$

where the CFL number is $C = \frac{c Δ t}{Δ x}$

$\tilde{ω} (k, c, Δ x, Δ t)$ is the numerical dispersion relation
$ω = k c$ is the exact dispersion relation
We can compare the two to investigate numerical accuracy and stability

Stability

A simpler approach is to insert for $u_{j}^{n} = A^{n} e^{i k x_{j}}$ directly in

$\frac{u_{j}^{n + 1} - 2 u_{j}^{n} + u_{j}^{n - 1}}{Δ t^{2}} = c^{2} \frac{u_{j + 1}^{n} - 2 u_{j}^{n} + u_{j - 1}^{n}}{Δ x^{2}}$

and solve for $A$ . We get

$\frac{(A^{n + 1} - 2 A^{n} + A^{n - 1}) e^{i k x_{j}}}{Δ t^{2}} = c^{2} A^{n} \frac{e^{i k (x_{j} + Δ x)} - 2 e^{i k x_{j}} + e^{i k (x_{j} - Δ x)}}{Δ x^{2}}$

Divide by $A^{n} e^{i k x_{j}}$ , multiply by $Δ t^{2}$ and use $C = c Δ t / Δ x$ to get

$A - 2 + A^{- 1} = C^{2} (e^{i k Δ x} - 2 + e^{- i k Δ x})$

continue on next slide

Stability

$A + A^{- 1} = 2 + C^{2} (e^{i k Δ x} - 2 + e^{- i k Δ x})$

Use $e^{i x} + e^{- i x} = 2 \cos x$ to obtain

$A + A^{- 1} = 2 + 2 C^{2} (\cos k Δ x - 1)$

This is a quadratic equation to solve for A. Using $β = 2 (1 + C^{2} (\cos (k Δ x) - 1))$ we get that

$A = \frac{β \pm \sqrt{β^{2} - 4}}{2}$

We see that $| A | = 1$ for any real numbers $- 2 \leq β \leq 2$ .

For all real numbers $- 2 \leq β \leq 2$

$| β \pm \sqrt{β^{2} - 4} | = 2$ since $| β \pm \sqrt{β^{2} - 4} | = | β + i \sqrt{4 - β^{2}} | = \sqrt{β^{2} + 4 - β^{2}} = 2$

For $| A | \leq 1$ and stability we need $- 2 \leq β \leq 2$ and thus

$- 2 \leq 2 (1 + C^{2} (\cos (k Δ x) - 1)) \leq 2$

Rearrange to get that

$- 2 \leq C^{2} (\cos (k Δ x) - 1) \leq 0$

Since $\cos (k Δ x)$ can at worst be $- 1$ we get that the positive real CFL number must be smaller than 1

$C \leq 1$

Hence (since $C = c Δ t / Δ x$ ) for stability we require that

$Δ t \leq \frac{Δ x}{c}$

Test Dirichlet solver using CFL=1.01 vs CFL=1.0

unm1[:] = sp.lambdify(x, u0.subs(t, 0))(xj)
un[:] = sp.lambdify(x, u0.subs(t, dt))(xj) 
plotdata = {0: unm1.copy()}
CFL = 1.01
for n in range(Nt):
  unp1[:] = 2*un - unm1 + CFL**2 * D2 @ un 
  unp1[0] = 0
  unp1[-1] = 0
  unm1[:] = un 
  un[:] = unp1
  if n % 10 == 0:
    plotdata[n] = unp1.copy()