Welcome to this bonus notebook that ties into the second module of "Practical Numerical Methods with Python".
In the first notebook of the second module, we were introduced to the numerical solutions of partial differential equations (PDEs). The equation used in that notebook was the 1-D linear convection equation, AKA the advection equation.
However, in this notebook we will explore the 1-D wave equation shown below:
\begin{equation}\frac{\partial^2 u}{\partial t^2}= c^2 \frac{\partial^2 u}{\partial x^2}. \end{equation}The wave equation is second order in both spatial and temporal dimensions.
If the initial displacement $\ f(x)$ and a velocity $\ g(x)$ are specified for the wave equation, then the solution of the system is: \begin{equation}\ u(x,t) = \frac{1}{2}[f(x-ct)+f(x+ct)]+\frac{1}{2c}\int{g(\tau){\partial \tau}}\end{equation} This solution is known as the d'Alembert solution for the 1-D wave equation.
For this notebook we will focus on the special case where the initial velocity is zero, $g(x)=0$ and the solution is simplified to the equation below: \begin{equation}\ u(x,t)=\frac{1}{2}[f(x-ct)+f(x+ct)]\end{equation}
A classical example of the 1-D wave equation is when simulating the motion of a guitar string after it is plucked.
As the figure above shows, the original triangular displacement splits into two waves traveling in opposite directions. One can also see that each of the travelling waves is half the height of the original function, which corresponds to d'Alembert's solution shown previously.
Now that we have established the solution for the 1-D wave equation, we will use this to simulate the longitudinal wave motion on a finite length bar.
Let's say we have a bar of finite length as shown below. The length is $10 \rm{m}$ in the $x$ direction. As can be seen, the system ends with a dashpot, also known as a damper; this dashpot is an important part of this problem because we will see how this will create a non-reflecting boundary condition.
In addition to the problem description, the initial displacement $f(x)$ is defined below: \begin{equation} u(x,0)=\begin{cases}2-2(x-5)^2 & \text{where } 4\leq x \leq 6,\\ 0 & \text{everywhere else in } (0, 10) \end{cases} \end{equation}
Since the 1-D wave equation is second order in time and space, we will discretize the equation using central difference for both dimensions. Below is the discretized equation:
Solving for the unknown in the discretized equation we get the following:
In this equation $C= c\frac{\Delta t}{\Delta x}$, which equates to the CFL number.
In addition to discretizing the wave equation, we also want to take a look at our boundary conditions. On the left-hand-side boundary, the bar is fixed, which would correlate with a Dirichlet boundary where the displacement is $u=0$.
The right-hand-side boundary is more tricky and that is where the idea of the nonreflecting boundary condition comes into play.
If we look at a simplified version of d'Alembert's solution, $f(x-ct)+f(x+ct),$ we can see that the first term represents the wave that is moving towards the right boundary while the latter term represents the wave that bounces away from the right boundary. In the bar we are analyzing, the dashpot at the right-hand-side means that the second return wave term would go to zero.
Differentiating the first term of D'Alembert's solution in terms of space and time would result in the following PDEs: \begin{equation}\frac{\partial u}{\partial x}= f^\prime,\frac{\partial u}{\partial t}= -cf^\prime\end{equation}
Rearranging the above equations gives the boundary condition we have to enforce on the right-hand side:
\begin{equation}\frac{\partial u}{\partial x}=-\frac{1}{c}\frac{\partial u}{\partial t}\end{equation}Discretizing our non-reflecting boundary:
and isolating for the end boundary:
The first thing we always do is to import the necessary libraries for this problem. In this case we will use bothy Numpy and Sympy. We will also use Matplotlib for plotting our displacement.
In [1]:
import numpy as np
import sympy as sp
from sympy.utilities.lambdify import lambdify
import matplotlib.pyplot as plt
%matplotlib inline
from matplotlib import rcParams
rcParams['font.family'] = 'serif'
rcParams['font.size'] = 16
We then create the initial condition for the bar.
In [2]:
# initial conditions and parameters
nx= 50
x= sp.symbols('x')
func= 2-2*(x-5)**2
init_func= lambdify((x), func)
x_array= np.linspace(0,10,nx+1)
u_initial= np.asarray([init_func(x0) for x0 in x_array])
Lets plot the initial function so we can start off on the right foot.
In [3]:
plt.plot(x_array, u_initial, color='#003366', ls='--', lw=3)
plt.xlim(0,10)
plt.ylim(0,3)
plt.xlabel('x')
plt.ylabel('Displacement')
Out[3]:
Now we can define a function to solve for the displacement.
In [4]:
def barsolver(u_init, T, C, nx):
''' Returns the displacement of the wave equation for a bar.
Parameters
----------
u_init: array of float
initial displacement of bar
T: integer
final time for calculation
C: integer
CFL number for stability
nx: integer
number of gridsteps
Returns
-------
u: array of float
final displacement of bar wave motion
'''
#initial parameters
c= 2
C2= C**2
dx= 10./nx
dt = C*dx/c
nt= int(round(T/dt))
#create arrays
u= np.zeros(nx+1) #array holding u for u[n+1]
u_1= np.zeros(nx+1) #array holding u for u[n]
u_2= np.zeros(nx+1) #array holding u for u[n-1]
u_1[4/dx:6/dx+1]= u_init[4/dx:6/dx+1]
# Loop for first time step
n = 0
for i in range(1, nx):
u[i] = u_1[i] + 0.5*C2*(u_1[i-1] - 2*u_1[i] + u_1[i+1])
# Enforce boundary conditions
u[0] = 0; u[-1] = -dx/(c*dt)*(u[i]-u_1[i])+u[-2]
# Switch variables before next step
u_2[:], u_1[:] = u_1, u
# Loops for subsequent time steps
for n in range(1, nt):
for i in range(1, nx):
u[i] = - u_2[i] + 2*u_1[i] + C2*(u_1[i-1] - 2*u_1[i] + u_1[i+1])
# Enforce boundary conditions
u[0] = 0; u[-1] = -dx/(c*dt)*(u[i]-u_1[i])+u[-2]
# Switch variables before next step
u_2[:], u_1[:] = u_1, u
return u
Lets run the solver for $t=0$ and plot the displacement.
In [5]:
u_ans= barsolver(u_initial,0,1,50)
In [6]:
plt.plot(x_array,u_ans)
plt.xlim(0,10)
plt.ylim(-3,3)
Out[6]:
As we expected, the displacement is equal to the initial parabola.
In [7]:
u_ans= barsolver(u_initial,2,1,50) #solving at t= 2
In [8]:
plt.plot(x_array,u_ans)
plt.xlim(0,10)
plt.ylim(-2,2)
Out[8]:
In [9]:
u_ans= barsolver(u_initial,4,1,50) #solving at t= 4
In [10]:
plt.plot(x_array,u_ans)
plt.xlim(0,10)
plt.ylim(-2,2)
Out[10]:
As we can see from the plots, the wave moving in the right direction has been absorbed while the left moving wave has flipped over and is now moving in the right direction.
Our code has expressed correctly the non-reflecting boundary!
What do you think will happen as more time passes?
You will need to plot the displacement for t= 6, 8, 10 and check if what you expected was correct.
Barba, Lorena A., et al. "MAE 6286 Practical Numerical Methods with Python," GW Open edX, The George Washingtion University, 2014. http://openedx.seas.gwu.edu/courses/GW/MAE6286/2014_fall/about .
Chapra, Steven C and Canale, Raymond P., "Numerical Methods for Engineers," McGraw-Hill Education, Inc., New York, 2010.
Everstine, Gordon C. "Analytical Solution of Partial Differential Equations," George Washington University, 2012. PDF
Fitzpatrick, Richard. "Introduction: Wave Equation," University of Texas, 2006. Wave Equation.
Strauss, Walter A., "Partial Differential Equations: An Introduction," Wiley, Inc., New York, 2008.
In [1]:
from IPython.core.display import HTML
css_file = '../../styles/numericalmoocstyle.css'
HTML(open(css_file, "r").read())
Out[1]: