Welcome to No Speeding Allowed, an individual contribution to Practical Numerical Methods with Python focusing on solution methods for hyperbolic PDEs by Chase Johnson. In this first module we will focus on an explicit finite difference method. In previous discussions, we've learned that explicit schemes have stability requirements but they are easier to understand and this is how we will be introducing second-order hyperbolic PDEs. In the following sections we will be building on the information we learn here to incorporate damping and then expand this to two dimensions.
We know that the general form of a second-order PDE is
$$a\frac{\partial^2 u}{\partial x^2}+b\frac{\partial^2 u}{\partial x y}+c\frac{\partial^2 u}{\partial y^2}+d\frac{\partial u}{\partial x}+e\frac{\partial u}{\partial y}+fu=g$$
If $b^2-4 a c >0$, the PDE is defined as hyperbolic and has wave-like solutions. This is the reason that hyperbolic PDEs are often referred to as "wave equations."
One of the most important aspects of wave equations that distinguish them from parabolic PDEs (Heat Equation) and elliptical PDEs (Laplace's Equation) is that information or disturbances are transmitted through space-time at a defined speed, hence the title of this lesson. I encourage you to read further on hyperbolic PDEs and Wikipedia has an excellent article.
We have already seen that looking at diffusion will yield parabolic PDEs, but where would we find natural phenomena that would yield hyperbolic PDEs?
The most well known second-order hyperbolic equation is the wave equation,
$$\frac{\partial^2 u}{\partial t^2} = c^2\frac{\partial^2 u}{\partial x^2}$$
and we will study variations of this equation in this section. Dr. Joel Feldman from the University of British Columbia has a very well written derivation of the wave equation that you can find here. You can also find derivations in References [Asmar, Brown].
Let's say that we have a string of length $L$ that is fixed at both ends and has a wave speed of $c$. We've already decided what our boundary conditions are, but what about our initial conditions? What options could we have?
The most common is to define an initial position and an initial velocity so our final problem statement is
$$\frac{\partial^2 u}{\partial t^2} = c^2\frac{\partial^2 u}{\partial x^2}$$
With boundary conditions: $u(0,t) = u(L,t) = 0$
And initial conditions: $u(x,0) = f(x)$ and $\frac{\partial u}{\partial t} = g(x)$
We have two second-order terms in our governing equation. We also have a first-order term in in our initial conditions, but we'll get to that later. We need to match our numerical scheme with the physics of our problem. We know that waves can travel in either direction, so we don't want to introduce a numerical bias that is only forward or backwards. We will start by discretizing our governing equation with central difference in both time and space.
Recall that the central difference approximation for a second-order derivative is:
$$\frac{d^2 u}{ds^2} \approx \frac{u_{i+1}-2u_i+u_{i-1}}{\Delta s^2}$$
Using this definition, let's rewrite our discretized governing equation and solve for $u_i^{n+1}$. We'll step through it together the first time.
Initially rewriting the derivative terms will give us:
$$\frac{u_i^{n+1}-2u_i^n+u_i^{n-1}}{\Delta t^2} = c^2\frac{u_{i+1}^n-2u_i^n+u_{i-1}^n}{\Delta x^2}$$
If we make the following substiution, $\lambda = \frac{c\Delta t}{\Delta x}$, and solve for $u_i^{n+1}$, we get:
$$u_i^{n+1} = 2u_i^n(1-\lambda^2)+\lambda^2(u_{i+1}^n+u_{i-1}^n)-u_i^{n-1}$$
Be sure that you can get the same result before moving on.
If you notice, our discretized governing equation will require values from the previous two time steps meaning that our method is not self starting. Just like Congress, we are going to have to give it a kickstart to get any meaningful work out of it. But unlike Congress, our method will continue to work until we tell it to stop!
Let's take a look at our $t=0$ array. This is where our initial position will be stored and we can easily accomplish this by simply forcing $u(x,0)=f(x)$.
But this only takes care of the first array, and our method will require not only $u(x,0)$ to be populated but also $u(x,1)$ to be populated before we can calculate $u(x,2)$.
Hmmm...what to do, what to do. Well, we still have an initial condition that we haven't used yet,
$$\frac{\partial u(x,0)}{\partial t}=g(x)$$
Using another finite difference, we can approximate
$$u(1,x)=u(0,x)+\Delta x g(x)$$
Now we have numerical solutions for our first and second time steps.
But...we have a problem. Did you notice that our approximation for the initial condition is only first-order accurate and our discretized scheme is second-order accurate? This means that our second time step will be less accurate and will feed our iterative loop propagating our error. We can do better than this, right?
We can derive a $\mathcal{O}(\Delta t^2)$ approximation for our first time step by considering the series expansion of $u(1,x)$, our second time step around our first time step, $u(0,x)$, which by definition would be a Maclaurin expansion.
$$u(1,x)=u(0,x)+\Delta t\frac{\partial}{\partial t} u(0,x)+\frac{\Delta t^2}{2}\frac{\partial^2}{\partial t^2} u(0,x)+\mathcal{O}(\Delta t^3)$$
Now we have our $\mathcal{O}(\Delta t^2)$ expression for $u(1,x)$, but we need to find expressions for the terms
$$\frac{\partial}{\partial t} u(0,x)\text{ and }\frac{\partial^2}{\partial t^2} u(0,x)$$
Using our governing equation we can show that
$$\frac{\partial^2}{\partial t^2}u(0,x) = c^2f"(x)$$
and by definition
$$\frac{\partial}{\partial t}u(0,x)=g(x)$$
Finally, we can get an expression for $u(1,x)$ that is $\mathcal{O}(\Delta t^2)$ accurate. Derive this yourself. $\mathit{Hint}$: Since we don't want to assume that $f"(x)$ necessarily exists, break it into another central difference approximation. You will find:
$$u(1,x_i) = (1-\lambda^2)f(x_i)+\frac{\lambda^2}{2}(f(x_{i+1})+f(x_{i-1})+\Delta t g(x)$$
All that is left is to program a loop to start populating the third time step based on the first two and to continue until our final time! There is one last bit to discuss before we get started. We are working with an explicit scheme, so of course Von Neumann stability analysis needs to be discussed. Dr. Peter Olver from the University of Minnesota has an excellent breakdown of both a graphical argument and von Neumann analysis that can be found here.
We will discuss the graphical representation here. Looking at the figure below from [Olver], we can see that time dependent values must remain inside the numerical domain of the mesh. If the wave speed is too high, we eventually would pull from values that have not been calculated yet. From this graphical representation we can see that
$$0 \leq c \leq \frac{\Delta x}{\Delta t} \text{ or } \Delta t \leq \frac{\Delta x}{c}$$
Now we are ready to begin!
First, let's try to model a guitar string plucked from an at rest position. This can be modeled with the following problem:
$$\frac{\partial^2 u}{\partial t^2}=64\frac{\partial^2 u}{\partial x^2}$$
$$u(0,t) = u(10,t) = 0$$
$$u(x,0) = 0$$
$$\frac{\partial u(x,0)}{\partial t} = 12cos(\frac{\pi(x-5)}{10})$$
In [10]:
# Importing our libraries
import numpy as np
import matplotlib.pyplot as plt
In [11]:
#Define Initial Parameters
L = 10. #Length of the string (given in problem statement)
nx = 100 #Number of steps in space
x = np.linspace(0,L,nx) #Useful vector for plotting
dx = L/(nx-1) #Size of steps in space
nt = 500 #Number of steps in time
c = 8.0 #Wave speed (given in problem statement)
dt = 0.5*dx/c #This comes from the CFL Condition
In [12]:
#Define Initial Conditions
yi = np.zeros(nx)
#Initial position profile (given in problem statement)
veli = 12*np.cos(np.pi*(x-5.)/L)
#Initial velocity profile (given in problem statement)
lam = c*dt/dx
In [13]:
#Plot Initial Conditions
plt.plot(x,yi,label='Initial Position')
plt.plot(x,veli,label='Initial Velocity')
plt.xlim(0,L)
plt.ylim(-0.5,12.5)
plt.legend(loc='best');
This is arguably one of the most important steps. It is very easy to get wrapped up in creating a function for our numerical model that we neglect the data that we are putting in. Before we get too deep into our problem, it is always a good idea to make sure that the intial conditions that we have created match our problem statement. Regardless of how good our scheme is, if our initial conditions aren't properly defined, our beautiful scheme won't work anyway. These look great, let's move on!
In [14]:
#Defining our Central Time/Central Space Scheme
def ctcs(nt,dt,ui,veli,lam):
'''Solves the wave equation with central-time, central-space scheme
Parameters:
--------------
nt: int
Number of time steps
dt: float
Size of time steps
ui: array of float
Initial position profile
veli: array of float
Initial velocity profile
lam: float
Wave constant defined as c*dt/dx
Returns:
--------------
y: array of float
Position profile of each time step stored in matrix. To call individual time step: y[i,:]
'''
y = np.zeros((nt,len(ui))) #Creates matrix to store values
y[0,:] = ui
y[1,1:-1] = (1-lam**2)*y[0,1:-1]+lam**2/2.*(y[0,2:]+y[0,:-2])+dt*veli[1:-1]
for i in range(1,nt-1):
y[i+1,1:-1] = 2.*y[i,1:-1]*(1-lam**2)+lam**2.*(y[i,2:]+y[i,:-2])-y[i-1,1:-1]
return y
In [15]:
#Let's run the program and plot a few solutions
y = ctcs(nt,dt,yi,veli,lam)
plt.figure(2)
plt.plot(x,y[0,:], label='Time 0')
plt.plot(x,y[4,:], label='Time 5')
plt.plot(x,y[9,:], label='Time 10')
plt.legend()
for i in range(nt):
plt.figure(3)
plt.plot(x,y[i,:])
plt.xlim(0,L)
That looks nice, but we have a better way to see all the solutions.
In [16]:
from matplotlib import animation
from JSAnimation.IPython_display import display_animation
def animate(data):
x = np.linspace(0,L,nx)
y = data
line.set_data(x,y)
return line,
fig = plt.figure();
ax = plt.axes(xlim=(0,L),ylim=(-6,6))
line, = ax.plot([],[],lw=2);
anim = animation.FuncAnimation(fig, animate, frames = y, interval = 50)
display_animation(anim,default_mode='reflect')
Out[16]: