Content is licensed under a Creative Commons Attribution 4.0 International License, code under MIT license $(c)$2014 Chase Johnson.

No Speeding Allowed

Solution Methods for 1-D Second-Order Hyperbolic PDEs

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.

Hyperbolic PDEs Defined

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.

Where do we see Hyperbolic PDEs

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

Defining the Problem

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

Explicit Finite Difference

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

CFL Condition for the Wave Equation

Now we are ready to begin!

The Vibrating String

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

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!

That looks nice, but we have a better way to see all the solutions.

That's pretty awesome. We can also model other wave phenomena with the same governing equation.

The Pulse Wave

Now let's try to model a traveling perturbation. We'll start by disturbing our string in the center of it's rest position and see how this information/energy is transmitted. This can be modeled with the following problem:
$$\frac{\partial^2 u}{\partial t^2}=9\frac{\partial^2 u}{\partial x^2}$$
$$u(0,t) = u(10,t) = 0$$
$$u(x,0) = 2 e^{-(x-5)^2}$$
$$\frac{\partial u(x,0)}{\partial t} = 0$$

What is really neat and unique about this wave equation is that information must travel at a finite speed. The endpoints of the string have no idea that there is a pulse until the information travels to the ends. This is unlike the diffusion equation that we were working with where a change in the boundaries was "felt" everywhere.

Dig Deeper...

In the next section you will be writing a code to model a more complete wave equation...one that models damping.
$$\frac{\partial^2 u}{\partial t^2}=c^2\frac{\partial^2 u}{\partial x^2}-\gamma\frac{\partial u}{\partial t}$$
How is this additional derivative term going to affect our technique? Make sure that you check out the second notebook in this series!

References:

• Ames, William F. "Numerical Methods for Partial Differential Equations." Third Edition. Academic Press, Inc. 1992.
• Asmar, Nakhle H. "Partial Differential Equations with Fourier Series and Boundary Value Problems." Second Edition. Pearson Education. 2005.
• Barba, Lorena A., et al. "MAE 6286 Practical Numerical Methods with Python." The George Washington University. http://openedx.seas.gwu.edu/courses/GW/MAE6286/2014_fall/about. 2014.
• Brown, James Ward and Ruel V. Churchill. "Fourier Series and Boundary Value Problems." Eighth Edition. McGraw-Hill. 2012.
• Faires, J. Douglas and Richard Burden. "Numerical Methods." Third Edition. Brooks/Cole. 2003.
• Feldman, Joel. "Derivation of the Wave Equation." http://www.math.ubc.ca/~feldman/m256/wave.pdf.
• Olver, Peter J. "Lecture Notes on Numerical Analysis." http://www.math.umn.edu/~olver/num.html.