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
# 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)
Failed to display Jupyter Widget of type HBox
.
If you're reading this message in the Jupyter Notebook or JupyterLab Notebook, it may mean that the widgets JavaScript is still loading. If this message persists, it likely means that the widgets JavaScript library is either not installed or not enabled. See the Jupyter Widgets Documentation for setup instructions.
If you're reading this message in another frontend (for example, a static rendering on GitHub or NBViewer), it may mean that your frontend doesn't currently support widgets.
Failed to display Jupyter Widget of type Figure
.
If you're reading this message in the Jupyter Notebook or JupyterLab Notebook, it may mean that the widgets JavaScript is still loading. If this message persists, it likely means that the widgets JavaScript library is either not installed or not enabled. See the Jupyter Widgets Documentation for setup instructions.
If you're reading this message in another frontend (for example, a static rendering on GitHub or NBViewer), it may mean that your frontend doesn't currently support widgets.