The 3D wave equation

Simulation of the three-dimensional wave equation using the finite difference method in Python.

The basics of the finite difference method

Simulation of standing waves by numerically solving the three-dimensional wave equation in Python.

Interference and diffraction of a wavefront at two circular holes.

After discussing the two-dimensional wave equation in an earlier article, lets now talk about its etension to three spatial dimensions. The principle is the same as with the two-dimensional wave equation, only now in three space dimensions.

This is the three-dimensional wave equation:

\begin{equation} \frac{\partial^2 u}{\partial t^2}(x,y,z,t) - \alpha^2 \left( \frac{\partial^2 u}{\partial x^2}(x,y,z,t) + \frac{\partial^2 u}{\partial y^2}(x,y,z,t) + \frac{\partial^2 u}{\partial z^2}(x,y,z,t) \right) = 0 \label{eq:pde1} \end{equation}

Those who have read the article on the 2-D wave equation already know that the actual purpose of this formulation is to scare students. One or the other prospective scientist may still suffer from a post-traumatic stress disorder from trying to solve the equation analytically because of the need to use the discrete Fourier transform in the process. If you are bold enough, you may read a description of the solution for only one space dimension in this link.

Jean Baptiste Joseph Fourier was undoubtedly a genius, for purely mathematically, his solution is an elegant and ingenious method, but that is not the point here, because this article is aimed at "non-geniuses" with a computer.

Equation \(\ref{eq:pde1}\) describes how a disturbance of the magnitude u moves in time and space in a medium in which waves propagate with the velocity α. To do this, it combines the second derivatives of the field u with respect to the spatial dimensions with the second derivative with respect to time. What exactly u is does not matter. It could stand for the air pressure, for example.

When physicists talk among themselves, they use a different formulation of the equation:

\begin{equation} u_{tt} - \alpha^2 \Delta u = 0 \label{eq:wave2d} \end{equation}

First of all it is noticeable that this is exactly the same equation as for the two-dimensional wave equation. The "triangle" in the equation stands for the Laplace operator, which is now three dimensional. It does not really matter whether 2D or 3D. In the three-dimensional case one must simply throw a couple more terms onto the problem to take the Z-Direction into account. For the solution we again need the second derivative of the field size with respect to time and the result of the application of the Laplace operator.

Since we are seeking a numerical solution, we need to discretize the problem. The simulation domain will be a three-dimensional grid. So we need three location dimensions and one time dimension. However, for the time dimension, we are only interested in the current time step and the two previous time steps. We need three fields that store the state of the field size at the three time points: \(u^{(0)}, u^{(1)}, u^{(2)}\).

The time index 0 stands for the future, the index 1 for the current state and the index 2 for the previous state of the field. The nodes of the grid have distance \(h\) from each other in the spatial dimension, and the time between two successive time steps is called \(t\). Each of the grid nodes represents an approximate value of the field \(u(x,y,z)\) at the location \(u(i*h, j*h, k*h)\).

We start with the determination of the time derivative \(u_{tt}\) at a grid node i,j at time index 1. For the numerical calculation of this derivative we use the 2nd order difference quotient. This requires, in addition to the current value of the field protery \(u^{(1)}_{(i,j,k)}\), a value from the future \(u^{(0)}_{(i,j,k)}\) and one from the past \(u^{(2)}_{(i,j,k)}\).

\begin{equation} u^{(1)}_{{tt}_{i,j,k}} = \frac{u^{(0)}_{i,j,k} - 2 u^{(1)}_{i,j,k} + u^{(2)}_{i,j,k}}{t^2} \label{eq:timediff} \end{equation}

Of course, we do not yet know the value from the future because this is what we want to compute. In addition to the second time derivative we need the Laplace operator for computing the 2nd order spatial derivative in x, y and z direction. The version of the three-dimensional, discrete Laplace operator \(\mathbf{D}^3_{i,j,k}\) shown here uses the values of the field property field in the 6 direct neighbor cells in x, y and z direction with single weighting and subtracts the value of the field in the node i,j,k with sixfold weighting:

\begin{equation} \quad \mathbf{D}^3_{i,j,k}= \left\{ \begin{bmatrix} 0 & 0 & 0 \\ 0 & 1 & 0 \\ 0 & 0 & 0 \end{bmatrix}, \begin{bmatrix} 0 & 1 & 0 \\ 1 & -6 & 1 \\ 0 & 1 & 0 \end{bmatrix}, \begin{bmatrix} 0 & 0 & 0 \\ 0 & 1 & 0 \\ 0 & 0 & 0 \end{bmatrix} \right\} \end{equation}

That looks complicated, but it's just an extension of the two-dimensional Laplace operator by one dimension, because the calculation of the second derivatives should now also take into account the two adjacent cells in the Z-direction. What is written down here in matrix form is the convolution kernel of the discrete three-dimensional Laplace operator. What it means is that the following equation for its determination at this point applies to the grid node i,j,k:

\begin{equation} \Delta u^{(1)}_{i,j,k} = \frac{u^{(1)}_{i-1,j,k}+u^{(1)}_{i+1,j,k}+u^{(1)}_{i,j-1,k}+u^{(1)}_{i,j+1,k} +u^{(1)}_{i,j,k-1}+u^{(1)}_{i,j,k+1} - 6 u^{(1)}_{i,j,k}}{h^2} \label{eq:laplace} \end{equation}

If we now insert the equations \(\ref{eq:timediff}\) and \(\ref{eq:laplace}\) into equation \(\ref{eq:wave2d}\), we obtain a discrete form of the three-dimensional wave equation:

\begin{equation} \frac{u^{(0)}_{i,j,k} - 2 u^{(1)}_{i,j,k} + u^{(2)}_{i,j,k}}{t^2} - \alpha^2 * \left( \frac{u^{(1)}_{i-1,j,k}+u^{(1)}_{i+1,j,k}+u^{(1)}_{i,j-1,k}+u^{(1)}_{i,j+1,k} +u^{(1)}_{i,j,k-1}+u^{(1)}_{i,j,k+1} - 6 u^{(1)}_{i,j,k}}{h^2} \right)= 0 \end{equation}

Now we rearrange this equation to solve for \(u^{(0)}_{(i,j,k)}\), since we are interested in the field value in the future:

\begin{equation} u^{(0)}_{i,j, k} = \left( \frac{t \alpha}{h}\right)^2 \left( u^{(1)}_{i-1,j,k}+u^{(1)}_{i+1,j,k}+u^{(1)}_{i,j-1,k}+u^{(1)}_{i,j+1,k} +u^{(1)}_{i,j,k-1}+u^{(1)}_{i,j,k+1} - 6 u^{(1)}_{i,j,k} \right) + 2 u^{(1)}_{i,j,k} - u^{(2)}_{i,j,k} \label{eq:findiff} \end{equation}

With this equation we can calculate every grid node in \(u^{(0)}\) from the two previous time steps \(u^{(1)}\) and \(u^{(2)}\). For this we iterate over all nodes of the field. At the beginning of the calculation, the values in the fields \(u^{(1)}\) and \(u^{(2)}\) are given as initial conditions. This is an explicit Euler method of second order.

More accurate spatial derivatives

The explicit Euler method of 2nd order is not particularly good and is usually not used for serious calculations. A simple way to increase the accuracy of the calculation is to improve the error term of the spatial derivatives used by the Laplace operator. This can be achieved by using a Laplace operator with a larger convolution kernel. Instead of a version with 3x3x3 entries (see equation \(\ref{eq:laplace}\)), one can use one with 5x5x5, 7x7x7 or even 9x9x9. These versions of the convolution kernel of the discrete Laplace operator can easily be derived from the two-dimensional versions of higher order. To demonstrate this, here is the Laplace operator of 4th order:

\(O(h4)\): \begin{equation} \quad \mathbf{D}^3_{i,j,k}= \frac{1}{12} \left\{ \begin{bmatrix}0 & 0 & 0 & 0 & 0\\0 & 0 & 0 & 0 & 0\\ 0 & 0 & -1 & 0 & 0\\0 & 0 & 0 & 0 & 0\\0 & 0 & 0 & 0 & 0\end{bmatrix}, \begin{bmatrix}0 & 0 & 0 & 0 & 0\\0 & 0 & 0 & 0 & 0\\ 0 & 0 & 16 & 0 & 0\\0 & 0 & 0 & 0 & 0\\0 & 0 & 0 & 0 & 0\end{bmatrix}, \begin{bmatrix}0 & 0 & -1 & 0 & 0\\0 & 0 & 16 & 0 & 0\\-1 & 16 & -90 & 16 & -1\\0 & 0 & 16 & 0 & 0\\0 & 0 & -1 & 0 & 0\end{bmatrix}, \begin{bmatrix}0 & 0 & 0 & 0 & 0\\0 & 0 & 0 & 0 & 0\\ 0 & 0 & 16 & 0 & 0\\0 & 0 & 0 & 0 & 0\\0 & 0 & 0 & 0 & 0\end{bmatrix}, \begin{bmatrix}0 & 0 & 0 & 0 & 0\\0 & 0 & 0 & 0 & 0\\ 0 & 0 & -1 & 0 & 0\\0 & 0 & 0 & 0 & 0\\0 & 0 & 0 & 0 & 0\end{bmatrix} \right\} \end{equation}

The extension to three dimensions consists of adding a Z-dimension with 4 additional matrices and entering the same values in the center of the matrices as in the x- and y-direction of the middle matrix. The value in the center of the 5x5x5 matrix must be lowered from -60 to -90, because the sum of the elements in the matrices added for the Z-direction is equal to 30. A version of the difference equation \(\ref{eq:findiff}\) with a spatial derivative term of error order 4 looks like this:

\begin{equation} \begin{aligned} u^{(0)}_{i,j,k} = \left( \frac{t \alpha}{h}\right)^2 \Bigl( & -u^{(1)}_{i-2,j,k}+16u^{(1)}_{i-1,j,k}+16u^{(1)}_{i+1,j,k}-u^{(1)}_{i+2,j,k} \\ & -u^{(1)}_{i,j-2,k}+16u^{(1)}_{i,j-1,k}+16u^{(1)}_{i,j+1,k}-u^{(1)}_{i,j+2,k} \\ & -u^{(1)}_{i,j,k-2}+16u^{(1)}_{i,j,k-1}+16u^{(1)}_{i,j,k+1}-u^{(1)}_{i,j,k+2} \\ & -90u^{(1)}_{i,j,k} \Bigr) + 2 u^{(1)}_{i,j,k} - u^{(2)}_{i,j,k} \label{eq:findiff16} \end{aligned} \end{equation}

The solution of the two-dimensional wave equation using spatial derivatives up to error order 6 is demonstrated in the file waves_3d_improved.py of the accompanying GitHub archive.

waves_3d_improved.py

Using spatial derivatives of higher error order leads to a significant reduction of numerical artifacts, but cannot completely suppress them. Accuracy is bought with higher computational effort but this also reduces the stability of the method. Therefore the maximum possible time step is slightly lower in comparison with the lower order schemes.

A general problem is that only the accuracy of the partial differences of the spatial derivatives has been improved. The time derivative is still calculated with an explicit second-order method. One might think: Lets use another high order explicit scheme for the time derivative too! Unfortunately this will not work as the resulting differce equation will become completely unstable and not produce any meaningfull results. Higher order time derivatives would require the use of so called implicit finite difference schemes which are beyond the scope of this article.

The model used here has a, in Z-direction, linearly increasing wave propagation velocity and a zone of constant low velocity in the center (e.g., a cavity in seismic exploration). The influence of the accuracy of the spatial derivatives is shown in the direct comparison of different implementations with O(h^2) (left), O(h^4) (middle) and O(h^6) (right). In particular, the simulation on the left shows significantly more numerical dispersion in the result.

Absorbing boundary conditions

The simulation, as described so far, has a decisive problem. All waves are reflected at the edge. If new wave peaks are placed again and again, after a short time nothing can be seen anymore, because the wave energy cannot leave the simulation area. What is the reason for that?

A 3x3x3 Laplace operator always needs a value in all directions next to the grid point to be calculated for the calculation of the Orstableitungen. For this reason, it cannot be applied to the grid points at the edge. We have therefore so far predetermined the edge values to 0. These are so-called Dirichlet boundary conditions. Therefore, no energy can be transferred outwards, which is why incoming waves are reflected back.

In order to simulate an open edge, we have to solve special differential equations on it. Such conditions are also called radiating, absorbing or Mur boundary conditions. [4], [5] For radiating boundary conditions, the following differential equations apply at the ends of the simulation area in the X direction:

\begin{align} \begin{alignedat}{3} \frac{\partial u}{\partial x} && \biggr\rvert_{x=0} &= -&&\frac{1}{c} \frac{\partial u}{\partial t} \label{eq:mur1} \nonumber\\ \frac{\partial u}{\partial x} && \biggr\rvert_{x=N} &= &&\frac{1}{c} \frac{\partial u}{\partial t} \label{eq:mur2} \end{alignedat} \end{align}

Analogous equations apply in the Y and Z directions. This results in the following difference equations for the boundary grid points x = 0, x = N, y = 0, y = N, z = 0 and z = N:

\begin{align} \begin{alignedat}{3} u^{(0)}_{0,j,k} &= u^{(1)}_{1,j,k} &&+ \frac{\kappa-1}{\kappa+1} (u^{(0)}_{1,j,k} &&- u^{(1)}_{0,j,k}) \\ u^{(0)}_{N,j,k} &= u^{(1)}_{N-1,j,k} &&+ \frac{\kappa-1}{\kappa+1} (u^{(0)}_{N-1,j,k} &&- u^{(1)}_{N,j,k}) \\ u^{(0)}_{i,0,k} &= u^{(1)}_{i,1,k} &&+ \frac{\kappa-1}{\kappa+1} (u^{(0)}_{i,1,k} &&- u^{(1)}_{i,0,k}) \\ u^{(0)}_{i,N,k} &= u^{(1)}_{i,N-1,k} &&+ \frac{\kappa-1}{\kappa+1} (u^{(0)}_{i,N-1,k} &&- u^{(1)}_{i,N,k}) \\ u^{(0)}_{i,j,0} &= u^{(1)}_{i,j,1} &&+ \frac{\kappa-1}{\kappa+1} (u^{(0)}_{i,j,1} &&- u^{(1)}_{i,k,0}) \\ u^{(0)}_{i,j,N} &= u^{(1)}_{i,j,N-1} &&+ \frac{\kappa-1}{\kappa+1} (u^{(0)}_{i,j,N-1} &&- u^{(1)}_{i,j,N}) \end{alignedat} \end{align}

where N stands for the index of the last grid node and \(\kappa\) is calculated as follows from wave propagation velocity, spatial step size and temporal step size:

\begin{equation} \kappa = \alpha \frac{t}{h} \end{equation}

These difference equations are applied to all grid points on the boundary. When using a larger discrete Laplace operator (5x5x5, 7x7x7 or 9x9x9), the boundary area must be chosen wider so that all grid node points that the Laplace operator cannot cover are included. The implementation of radiating boundary conditions according to Mur is demonstrated in the file waves_3d_improved.py of the accompanying GitHub archive.

waves_3d_improved.py

Python Source Code

The short version of the source code shown here uses spatial derivatives with an error term of second order and 1st order absorbing boundary conditions according to Mur. It is based on that of the two-dimensional verison. The mayavi library is used for 3d rendering.

waves_3d.py

import numpy as np
import math

from mayavi import mlab
from tvtk.util.ctf import *

hs = ts = 1        # time and space step width
dimx = dimy = 100  # dimensions of the simulation domain 
dimz = 60           

u = np.zeros((3, dimx, dimy, dimz))   # three state arrays of the field variable

# Initial conditions
vel = np.full((dimx, dimy, dimz), 0.3)
vel[0:dimx, 0:dimy, 4:6] = 0.0

for i in range(dimx):
    for j in range(dimy):
        for k in range(4, 6):
            distance1 = np.linalg.norm(np.array([i, j, k]) - np.array([40, 50, 5]))
            distance2 = np.linalg.norm(np.array([i, j, k]) - np.array([60, 50, 5]))
            if distance1 <= 5 or distance2 <= 5:
                vel[i, j, k] = 0.3

tau = ( (vel*ts) / hs )**2
k = ts * vel / hs  

def update(u : any):
    u[2] = u[1]
    u[1] = u[0]

    # the 3D Wave equation
    u[0, 1:dimx-1, 1:dimy-1, 1:dimz-1] = tau[1:dimx-1, 1:dimy-1, 1:dimz-1] \
        * (   u[1, 0:dimx-2, 1:dimy-1, 1:dimz-1] +     u[1, 1:dimx-1, 0:dimy-2, 1:dimz-1]
            + u[1, 1:dimx-1, 1:dimy-1, 0:dimz-2] - 6 * u[1, 1:dimx-1, 1:dimy-1, 1:dimz-1]
            + u[1, 2:dimx  , 1:dimy-1, 1:dimz-1] +     u[1, 1:dimx-1, 2:dimy,   1:dimz-1]
            + u[1, 1:dimx-1, 1:dimy-1, 2:dimz] ) + 2 * u[1, 1:dimx-1, 1:dimy-1, 1:dimz-1] - u[2, 1:dimx-1, 1:dimy-1, 1:dimz-1]

    # Absorbing boundary conditions for all 6 sides of the simulation domain
    c = dimx - 1
    u[0, dimx-2:c, 1:dimy-1, 1:dimz-1] = u[1, dimx-3:c-1, 1:dimy-1, 1:dimz-1] \
        + (k[dimx-2:c, 1:dimy-1, 1:dimz-1]-1)/(k[dimx-2:c, 1:dimy-1, 1:dimz-1]+1) \
        * (u[0, dimx-3:c-1, 1:dimy-1, 1:dimz-1] - u[1, dimx-2:c,1:dimy-1, 1:dimz-1])
    c = 0
    u[0, c:1, 1:dimy-1, 1:dimz-1] = u[1, 1:2, 1:dimy-1, 1:dimz-1] \
        + (k[c:1, 1:dimy-1, 1:dimz-1]-1)/(k[c:1, 1:dimy-1, 1:dimz-1]+1) \
        * (u[0, 1:2,1:dimy-1, 1:dimz-1] - u[1, c:1,1:dimy-1, 1:dimz-1])
    r = dimy-1
    u[0, 1:dimx-1, dimy-2:r, 1:dimz-1] = u[1, 1:dimx-1, dimy-3:r-1, 1:dimz-1] \
        + (k[1:dimx-1, dimy-2:r, 1:dimz-1]-1)/(k[1:dimx-1, dimy-2:r, 1:dimz-1]+1) \
        * (u[0, 1:dimx-1, dimy-3:r-1, 1:dimz-1] - u[1, 1:dimx-1, dimy-2:r, 1:dimz-1])
    r = 0
    u[0, 1:dimx-1, r:1, 1:dimz-1] = u[1, 1:dimx-1, 1:2, 1:dimz-1] \
        + (k[1:dimx-1, r:1, 1:dimz-1]-1)/(k[1:dimx-1, r:1, 1:dimz-1]+1) \
        * (u[0, 1:dimx-1, 1:2, 1:dimz-1] - u[1, 1:dimx-1, r:1, 1:dimz-1])
    d = dimz-1
    u[0, 1:dimx-1, 1:dimy-1, dimz-2:d] = u[1, 1:dimx-1, 1:dimy-1, dimz-3:d-1] \
        + (k[1:dimx-1, 1:dimy-1, dimz-2: d]-1)/(k[1:dimx-1, 1:dimy-1, dimz-2:d]+1) \
        * (u[0, 1:dimx-1, 1:dimy-1, dimz-3:d-1] - u[1, 1:dimx-1, 1:dimy-1, dimz-2:d])
    d = 0
    u[0, 1:dimx-1, 1:dimy-1, d:1] = u[1, 1:dimx-1, 1:dimy-1, d+1:2] \
        + (k[1:dimx-1, 1:dimy-1, d:1]-1)/(k[1:dimx-1, 1:dimy-1, d:1]+1) \
        * (u[0, 1:dimx-1, 1:dimy-1, d+1:2] - u[1, 1:dimx-1, 1:dimy-1, d:1])

@mlab.animate(delay=20)
def update_loop():
    fig = mlab.figure(bgcolor=(0, 0, 0))
    src = mlab.pipeline.scalar_field(u[0])
    volume = mlab.pipeline.volume(src)
    tick = 0
    while True:
        tick += 0.1
        mlab.view(azimuth=4*tick, elevation=80, distance=200, focalpoint=(int(dimx/2), int(dimy/2), 20))

        u[0:2, 1:dimx-1, 1:dimy-1, 1:3] += 2 * math.sin(tick*1.5)
        update(u)

        src.mlab_source.scalars = u[0]
        absmax = 1 

        # define opacity function
        otf = PiecewiseFunction()
        for val, opacity in [(absmax, 1), (absmax * 0.2, 0), (-absmax, 1), (-absmax * 0.2, 0)]:
            otf.add_point(val, opacity)
        volume._volume_property.set_scalar_opacity(otf) 

        # define color function (red-black-green)
        ctf = ColorTransferFunction()
        for p in [(absmax, 0, 1, 0), (0, 0, 0, 0), (-absmax, 1, 0, 0)]:
            ctf.add_rgb_point(*p)
        volume._volume_property.set_color(ctf)

        yield

if __name__ == "__main__":
    animate = update_loop()
    mlab.show()

References

J. Douglas Faires, Richard L. Burden: "Numerische Methoden: Näherungsverfahren und ihre praktische Anwendung." Spektrum Akademischer Verlag; ISBN 3-86025-332-8; 1995
Hans Stefani, Gerhard Kluge: "Theoretische Mechanik" Spektrum Akademischer Verlag; ISBN 3-86025-284-4; 1995
Bengt Fornberg: "Generation of finite difference formulas on arbitrarily spaced grids." In: Mathematics of Computation. Band 51, Nr. 184, 1988, Seite 702
Gerrit Mur: "Absorbing Boundary Conditions for the Finite-Difference Approximation of the Time-Domain Electromagnetic-Field Equations" IEEE Transactions on Electromagnetic Compatibility: Vol EMC-23, November 1981; Page 377-382
unbekannter Autor: "Finite Difference Methods" online; unbekannter Herausgeber