The wave equation is a hyperbolic partial differential equation (PDE) which describes the displacement \(y(x,t)\) as a function of position and time. The 1D linear wave equation is

\[\frac{\partial^2 y(x, t)}{\partial x^2} = \frac{1}{c^2} \frac{\partial^2 y(x, t)}{\partial t^2}\]

whose solutions are traveling waves with wave velocity \(c\), for example, the waves that are generated by plucking a string:

The wave equation in 2D can be derived in the same fashion as for the 1D case and reads

\[\left(\frac{\partial^2}{\partial x^2} + \frac{\partial^2}{\partial y^2}\right)u(x, y, t) - \frac{1}{c^2} \frac{\partial^2 u}{\partial t^2} = 0\]

or in general (3D)

\[\boldsymbol{\nabla}^2 u(\mathbf{x}, t) - \frac{1}{c^2} \frac{\partial^2}{\partial t^2}u(\mathbf{x}, t) = 0.\]

Leap frog algorithm for solving the wave equation

We can then transform the PDE into a finite difference equation on a \(\Delta x, \Delta t\) lattice. The difference equation can be solved with a time stepping scheme where we start from the initial values and solve the spatial component for increasing times \(t = \Delta t, 2\Delta t, 3\Delta t, \dots\) using an explicit method such as the Leapfrog algorithm, similar to the solution for the heat equation.

For example, in the 2D case, the finite difference approximation for the Laplace operator is

\[\begin{split} \left(\frac{\partial^2}{\partial x^2} + \frac{\partial^2}{\partial y^2}\right)u(x, y, t) &\approx \frac{u(x+\Delta x, y, t) - 2u(x, y, t) + u(x-\Delta x, y, t)}{\Delta x^2} \\ &+ \frac{u(x, y + \Delta y, t) - 2u(x, y, t) + u(x, y-\Delta y, t)}{\Delta y^2} \end{split}\]

Setting \(\Delta x = \Delta y = \delta\) (same spacing in \(x\) and \(y\)) and \(u_{ijn} := u(i\delta, j\delta, n\Delta t)\) we obtain the discretized 2D wave equation

\[u_{i,j,n+1} = 2(1 - 2\beta^2)u_{i,j,n} - u_{i,j,n-1} + \beta^2(u_{i+1,j,n} + u_{i-1,j,n} + u_{i,j+1,n} + u_{i,j-1,n})\\ \beta := \frac{c}{\delta / \Delta t}\]

where the unknown (future) term is on the left hand site and can be computed from all known (present and past) terms on the right hand side.

Von Neumann stability analysis

Not all combinations of \(\Delta t\) and \(\Delta x\) lead to stable solutions. Using von Neumann stability analysis one can determine analytically the relationship between the discretizations in space and time that lead to stable solutions. The analysis is based on determining the stable eigenmodes of the finite difference equation, i.e., we determine the conditions under which no modes grow in time. This then implies that solutions, which can be written as linear superpositions of these eigenmodes, will also not grow in time and are therefore stable.

For the wave equation the Courant condition follows

\[c \le \frac{\Delta x}{\Delta t}\]

(1D case) or when written as

\[\Delta t \le \frac{\Delta x}{c}\]

it means that the time step \(\Delta t\) (for a given \(\Delta x\)) must be smaller than the time that the wave takes to travel one grid step.

In the 2D case, the Courant condition is \(c \le \frac{1}{\sqrt{2}} \frac{\delta}{\Delta t}\).

Class material

Additional resources

  • Computational Physics, Ch 21
  • Computational Modeling, Ch 6.5
  • Numerical Recipes in C, WH Press, SA Teukolsky, WT Vetterling, BP Flannery. 2nd ed, 2002. Cambridge University Press. Chapter 19.

Footnotes

  1. As usual, git pull the resources repository to get a local copy of the notebook. Then copy the notebook into your work directory in order to complete the exercises. ↩

  2. Notebook will be posted after class; in the mean time look at the student notebook. ↩