Welcome back to our next step in understanding hyperbolic PDEs, in particular the wave equation. We are going to increase the difficulty by altering our governing equation.
Recall that our previous governing equation was:
$$\frac{\partial^2 u}{\partial t^2}=c^2\frac{\partial^2 u}{\partial x^2}$$
For this assignment we are going to incorporate a damping term:
$$\frac{\partial^2 u}{\partial t^2}=c^2\frac{\partial^2 u}{\partial x^2}-\gamma\frac{\partial u}{\partial t}$$
Where $c$ is still the wave speed and now $\gamma$ is the damping coefficient. You can see that the larger $\gamma$ is, the faster the motion is going to be diminished.
Your challenge is to develop a code that can solve the following problem:
Let's try to model a disturbance in a resistive media. This will be somewhat like dropping a pebble into water, but in 1-D. We can also start our time (since starting point is arbitrary as long as we are consistent) right at the maximum point of the deflection so our initial velocity will be zero. With that in mind, here is your problem:
$$\frac{\partial^2 u}{\partial t^2}=64\frac{\partial^2 u}{\partial x^2}-5.5\frac{\partial u}{\partial t}$$
$$u(0,t)=u(10,t)=0$$
$$u(x,0)=-12e^{-4(x-5)^2}$$
$$\frac{\partial u(x,0)}{\partial t}=0$$
Be sure to use the same method to get our $\mathcal{O}(\Delta t^2)$ accuracy for $u(x,1)$. It will be slightly different than last time due to the new damping term. Use the second-order central difference for both second-order derivatives and foward difference in time for the damping term.
I've included the solution equations that I derived in a similar fashion to the previous notebook. The level of difficulty is slightly higher because now there is a first order time derivative to deal with on the right hand side. The student will still need to pre-fill the zeroth and first time step as we walked through before. The student will need to realize that $\frac{\partial u[0,x]}{\partial t}$ is really just $g(x)$ given in the problem statement.
$$u_i^1=u_i^0(1-\lambda^2)+\frac{\lambda^2}{2}(u_{i+1}^0+u_{i-1}^0)+\Delta t g(x)(1-\frac{\Delta t \gamma}{2})$$
$$u_i^{n+1}=\frac{2u_i^n(1-\lambda^2-\frac{\gamma \Delta t}{2})+\lambda^2(u_{i+1}^n+u_{i-1}^n)-u_i^{n-1}}{1+\gamma \Delta t}$$
In [3]:
import numpy as np
import matplotlib.pyplot as plt
#Define Initial Parameters
u_max = 2.
L = 10.
nx = 100
x = np.linspace(0,L,nx)
dx = L/(nx-1)
nt = 100
c = 8.0
dt = 0.5*dx/c
gamma = 5.5
In [4]:
yi = -12.*np.exp(-4.*(x-5.)**2)
veli = np.zeros(nx)
lam = c*dt/dx
#Plot Initial Conditions
plt.plot(x,yi)
plt.plot(x,veli)
plt.xlim(0,L)
plt.ylim(-12,0.1)
Out[4]:
In [11]:
def ctcs_damp(nt,dt,ui,veli,lam,gamma):
'''Solves the wave equation with central-time, central-space scheme with damping
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
gamma: float
Damping coefficient
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(yi)))
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]*(1-dt*gamma/2.)
for i in range(1,nt-1):
y[i+1,1:-1] = (2.*y[i,1:-1]*(1-lam**2-(gamma*dt/2.))+lam**2.*(y[i,2:]+y[i,:-2])-y[i-1,1:-1])/(1.+gamma*dt)
return y
Notice that the student can't easily just cut and paste from the last lesson. They need to recalculate the solution equations above by hand. Then it is not entirely difficult for them to adapt their previous ctcs code to accomodate the new equations. Here is the final animation.
In [13]:
y = ctcs_damp(nt,dt,yi,veli,lam,gamma)
In [14]:
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=(-15,10))
line, = ax.plot([],[],lw=2);
anim = animation.FuncAnimation(fig, animate, frames = y, interval = 50)
display_animation(anim,default_mode='loop')
Out[14]:
This is a very good exercise and is a logical step up from the previous walkthrough. The level of difficulty is moderate, but is doable. To assess the student the amplitude at a given location for a given time step could be used as a grading criterion.
In [15]:
from IPython.core.display import HTML
css_file = '../../styles/numericalmoocstyle.css'
HTML(open(css_file, "r").read())
Out[15]:
In [ ]: