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)$.

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). $$

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)
dt=1.0
c=1.0
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)
y[:N]=numpy.exp(-x**2)
Y=[]
for _ in range(Nt):
Y.append(y)
z=numpy.dot(B,y)
y = linalg.solve(A, z)
Y.append(y)
pyplot.figure(1)
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,
value=0,
min=0,
max=Nt,
step=1,
description="Press play",
disabled=False
)
ipywidgets.jslink((play, 'value'), (slider, 'value'))
display(ipywidgets.HBox([play, slider]))
display(plt)
```

In [ ]:

```
```

In [ ]:

```
```