The Wave Equation

The wave equation is one of the canonical partial differential equationss with which numericists deal, and an example of a hyperbolic PDE. As you probably know, the equation for a parabola is

$$ x^2 = a^2 y^2+b. $$

Similar to this, the standard wave equation in one dimension is

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

This is a model for a system radiating information, where the speed of the propagation is $c$, whichand has units of m/s.

When $c$ is constant, the equation accepts leftward and rightward travelling wave solutions, $f(x+ct)$ and $g(x-ct)$.

Initial Conditions

The second order ODE can be rewritten as a pair of coupled first order ODEs, $$ \frac{\partial v}{\partial t} = c^2\frac{\partial ^2 u}{\partial t}, $$ $$ \frac{\partial u}{\partial t} = v. $$ where the second equation is a rewritten version of the identity, $v=\frac{\partial u}{\partial t}$. From this, it can be seen that we need to provide initial conditions for both $u$ and $v$, i.e $$ u(x,0)= U(x) $$ $$\frac{\partial u}{\partial t} (x,0)= V(x). $$

Boundary Conditions

The boundary conditions must supply enough information to control both leftward and rightward travelling waves, so that we need a boundary condition on every boundary.

Lets build a finite difference model for a wave equation

In [44]:
# import relevant modules

import numpy
from scipy import linalg
from bqplot import pyplot
from IPython.display import display
import ipywidgets

# set up our discretization

N = 101
Nt= 100
dx = 10.0/(N-1)
x = numpy.linspace(-10.0, 10.0, N)

# useful matrices. A handles the new step, B the old.

A = numpy.zeros([2*N, 2*N])
B = numpy.zeros([2*N, 2*N])

for j in range(N):
    # Dirichlet boundary conditions
    A[j,j] = 1.0
    A[N-1,N-1] = 1.0
    A[N+j,N+j] = 1.0
    A[2*N-1,2*N-1] = 1.0
for i in range(1,N-1):
    A[i, i] = 1.0
    A[i, N+i-1] =-dt*c*1.0/dx**2/2.0
    A[i, N+i] = +dt*c*2.0/dx**2/2.0
    A[i, N+i+1] = -dt*c*1.0/dx**2/2.0
    A[N+i, i] = -dt/2.0
    A[N+i, N+i] = 1.0
    B[i, i] = 1.0
    B[i, N+i-1] = dt*c*1.0/dx**2/2.0
    B[i, N+i] = -dt*c*2.0/dx**2/2.0
    B[i, N+i+1] = dt*c*1.0/dx**2/2.0
    B[N+i, i] = dt/2.0
    B[N+i, N+i] = 1.0

# initial conditions
y = numpy.zeros(2*N)

for _ in range(Nt):
    y = linalg.solve(A, z)


xsc = pyplot.LinearScale(min=-10,max=10)
ysc = pyplot.LinearScale(min=-1,max=1)
axy = pyplot.Axis(label='u', scale=ysc, orientation='vertical', side='left', grid_lines='solid')
axx = pyplot.Axis(label='x', scale=xsc, grid_lines='solid')
lines = pyplot.Lines(x=x, y=Y[0][:N], scales = {'x': xsc, 'y': ysc})
plt = pyplot.Figure(layout=ipywidgets.Layout(width='auto'), marks=[lines],axes=[axx,axy])
def g(t):
    lines.y = Y[int(t/dt)][:N]
w = ipywidgets.interactive(g, t=(0.,1.0*Nt,1));
slider = w.children[0]
output = w.children[-1]

play = ipywidgets.Play(
#     interval=10,
    description="Press play",

ipywidgets.jslink((play, 'value'), (slider, 'value'))

display(ipywidgets.HBox([play, slider]))
In [ ]:
In [ ]: