*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.

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

$$\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].

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

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

$$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?

$$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!

$$\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');

```
```

```
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]:
```