From COMPUTATIONAL PHYSICS, 3rd Ed, 2015 RH Landau, MJ Paez, and CC Bordeianu (deceased) Copyrights: Wiley-VCH, Berlin; and Wiley & Sons, New York R Landau, Oregon State Unv, MJ Paez, Univ Antioquia, C Bordeianu, Univ Bucharest, 2015. Support by National Science Foundation. |
21 Wave Equations I: Strings & Membranes
21.1 A Vibrating String
21.2 The Hyperbolic Wave Equation (Theory)
21.2.1 Solution via Normal-Mode Expansion
21.2.2 Algorithm: Time-Stepping
21.2.3 Wave Equation Implementation
21.2.4 Assessment, Exploration
21.3 Strings with Friction (Extension)
21.4 Strings with Variable Tension & Density
21.4.1 Waves on Catenary
21.4.2 Derivation of Catenary Shape
21.4.3 Catenary & Frictional Wave Exercises
21.5 Vibrating Membrane (2-D Waves)
21.6 Analytical Solution
21.7 Numerical Solution for 2-D Waves
In this chapter and the next we explore the numerical solution of several PDE’s that yield waves as solutions. If you have skipped the discussion of the heat equation in Chapter 20, PDE Review & Heat Flow via Time Stepping,* then this chapter will be the first example of how initial conditions are propagated forward in time with a time-stepping or leapfrog algorithm. An important aspect of this chapter is its demonstration that once you have a working algorithm for solving a wave equation, you can include considerably more physics than is possible with analytic treatments. First we deal with a number of aspects of 1-D waves on a string, and then with 2-D waves on a membrane. In the next chapter we look at quantum wave packets and E&M waves. In Chapter 24 we look at shock and solitary waves.*
** This Chapter’s Lecture, Slide Web Links, Applets & Animations**
All Lectures |
| Lecture (Flash)| Slides | Sections|Lecture (Flash)| Slides | Sections|
|- - -|:- - -:|:- - -:|- - -|:- - -:|:- - -:|
|Realistic String Waves|pdf| 18.1| Catenary Waves|pdf|18.3|
|Quantum WavePackets| pdf|18.5 | EM Waves (FDTD)|pdf|22.5|
|Cat with Friction Applet|-|21.4|Waves on a String Applet|-| 18.1-18.2 |
|String normal mode|-|18.1-18.2 |Two peaks on a String Applet|-| 18.1-18.2|
|Catenary Wave Movie|-|18.4|Wavepacket-Wavepacket Scattering Movies|-|18.5-18.7 |
|Wavepacket Slit Movie|-| 18.5-18.7| Square Well Applet|-|18.6 |
|Harmonic Oscillator|-|18.7|Two Slit Interference|-|18.8|
Problem: Recall the demonstration from elementary physics in which a string tied down at both ends is plucked “gently” at one location and a pulse is observed to travel along the string. Likewise, if the string has one end free and you shake it just right, a standing-wave pattern is set up in which the nodes remain in place and the antinodes move just up and down. Your problem is to develop an accurate model for wave propagation on a string, and to see if you can set up traveling- and standing-wave patterns.[Note: Some similar but independent studies can also be found in [Rawitscher et al.(96)].]
Consider a string of length
Figure 21.1 A: A stretched string of length
To obtain a simple linear equation of motion (nonlinear wave equations are
discussed in Chapters 24, Shock Waves & Solitons, and 25,
Fluid Dynamics) we assume that the string’s relative
displacement
$$\begin{align} \tag*{21.1} \sum F_y &= \rho \Delta x \frac{\partial^2 y}{\partial t^2}\ & = T\sin\theta(x+\Delta x) - T\sin\theta(x)\ & = T \left. \frac{\partial y}{\partial x}\right|_{x+\Delta x}
- T \left.\frac{\partial y}{\partial x}\right|_x \simeq T \frac{\partial^2 y}{\partial x^2}, \tag*{21.3} \ \Rightarrow \quad \frac{\partial^2 y(x,t)}{\partial x^2} & = \frac{1}{c^2} \frac{\partial^2 y(x,t)}{\partial t^2}, \quad c = \sqrt{\frac {T}{\rho}}, \tag*{21.4} \end{align}$$
where we have assumed
that
The initial condition for our problem is that the string is plucked
gently and released. We assume that the “pluck” places the string in a
triangular shape with the center of triangle
Because (21.4) is second-order in time, a second initial condition (beyond initial displacement) is needed to determine the solution. We interpret the “gentleness” of the pluck to mean the string is released from rest:
The boundary conditions have both ends of the string tied down for all times:
The analytic solution to (21.4) is obtained via the familiar separation-of-variables technique. We assume that the solution is the product of a function of space and a function of time:
We substitute (21.8) into (21.4), divide by
The angular frequency
The time solution is
where the frequency of this $n$th normal mode is also fixed. In fact,
it is the single frequency of oscillation that defines a normal mode.
The initial condition (21.5) of zero velocity,
Figure 21.2 The solutions of the wave equation for four earlier space-time points are used to obtain the solution at the present time. The boundary and initial conditions are indicated by the white-centered dots.
Because the wave equation (21.4) is linear in
(Yet we will lose linear superposition once we include nonlinear terms in the
wave equation.) The Fourier coefficient
We multiply both sides by
You will be asked to compare the Fourier series (21.14) to our numerical solution. While it is in the nature of the approximation that the precision of the numerical solution depends on the choice of step sizes, it is also revealing to realize that the precision of the so-called analytic solution depends on summing an infinite number of terms, which can be carried out only approximately.
As with Laplace’s equation and the heat equation, we look for a solution
In
contrast to Laplace’s equation where the grid was in two space dimensions, the
grid in Figure 21.2 is in both space and time. That being the case, moving across
a row corresponds to increasing
As with the Laplace equation, we use the central-difference approximation to discretize the wave equation into a difference equation. First we express the second derivatives in terms of finite differences:
Substituting (21.19) in the wave equation (21.4) yields the difference equation
Notice that this equation contains three time values: $j+1 = $ the future, $j = $ the present, and $j - 1 = $ the past. Consequently, we rearrange it into a form that permits us to predict the future solution from the present and past solutions:
Here
Figure 21.2 The solutions of the wave equation for four earlier space-time points are used to obtain the solution at the present time. The boundary and initial conditions are indicated by the white-centered dots.
As you have seen in our discussion of the heat equation, a leapfrog method is quite different from a relaxation technique. We start with the solution along the topmost row and then move down one step at a time. If we write the solution for present times to a file, then we need to store only three time values on the computer, which saves memory. In fact, because the time steps must be quite small to obtain high precision, you may want to store the solution only for every fifth or tenth time.
Initializing the recurrence relation is a bit tricky because it requires displacements from two earlier times, whereas the initial conditions are for only one time. Nonetheless, the rest condition (21.5) when combined with the central-difference approximation lets us extrapolate to negative time:
Here we take the initial time as
Equation (21.23) uses the solution throughout all space at the initial
time
As is also true with the heat equation, the success of the numerical method
depends on the relative sizes of the time and space steps. If we apply a von
Neumann stability analysis to this problem by substituting
Equation (21.24) means that the solution gets better with smaller time steps
but gets worse for smaller space steps (unless you simultaneously make the time
step smaller). Having different sensitivities to the time and space steps may
appear surprising because the wave equation (21.4) is symmetric in
Exercise: Figure out a procedure for solving for the wave equation for all times in just one step. Estimate how much memory would be required.
Figure 21.3 The vertical displacement as a function of position
Exercise: Try to figure out a procedure for solving for the wave motion with a relaxation technique. What would you take as your initial guess, and how would you know when the procedure has converged?
#### EqStringAnimate.py, Notebook Version, Vibrating string + MatPlotLib
from numpy import *
import numpy as np
import matplotlib.pyplot as plt
#import scipy.integrate as integrate
import matplotlib.animation as animation
# Parameters
rho = 0.01 # string density
ten = 40. # string tension
c = sqrt(ten/rho) # Propagation speed
c1 = c # CFL criterium
ratio = c*c/(c1*c1)
# Initialization
xi = np.zeros( (101, 3), float) # 101 x's & 3 t's
k=range(0,101)
def init():
for i in range(0, 81):
xi[i, 0] = 0.00125*i # Initial condition: string plucked,shape
for i in range (81, 101): # first part of string
xi[i, 0] = 0.1 - 0.005*(i - 80) # second part of string
init() # plot string initial position
fig=plt.figure() # figure to plot (a changing line)
# select axis; 111: only one plot, x,y, scales given
ax = fig.add_subplot(111, autoscale_on=False, xlim=(0, 101), ylim=(-0.15, 0.15))
ax.grid() # plot a grid
plt.title("Vibrating String")
line, = ax.plot(k, xi[k,0], lw=2) # x axis, y values, linewidth=2
# Later time steps
for i in range(1, 100): # use algorithm
xi[i, 1] = xi[i, 0] + 0.5*ratio*(xi[i + 1, 0] + xi[i - 1, 0] - 2*xi[i, 0])
def animate(num): #num: dummy, algorithm, will plot (x, xi)
for i in range(1, 100):
xi[i,2] = 2.*xi[i,1]-xi[i,0]+ratio*(xi[i+1,1]+xi[i-1,1]-2*xi[i,1])
line.set_data(k,xi[k,2]) # data to plot ,x,y
for m in range (0,101): # part of algorithm
xi[m, 0] = xi[m, 1] # recycle array
xi[m, 1] = xi[m, 2]
return line,
# next: animation(figure, function,dummy argument: 1
ani = animation.FuncAnimation(fig, animate,1)
plt.show()
print("finished")
The program EqStringAnimate.py
in Listing 21.1 solves the wave
equation for a string of length
Listing 21.1 EqStringAnimate.py solves the wave equation via time
stepping for a string of length
Figure 21.4 The vertical displacement as a function of position and time of a string initially plucked
simultaneously at two points, as shown by arrows. Note that each
initial peak breaks up into waves traveling to the right and to
the left. The traveling waves invert on reflection from the
fixed ends. As a consequence of these inversions, the
-
Solve the wave equation and make a surface plot of displacement versus time and position.
-
Explore a number of space and time step combinations. In particular, try steps that satisfy and that do not satisfy the Courant condition (21.24). Does your exploration confirm the stability condition?
-
Compare the analytic and numeric solutions, summing at least 200 terms in the analytic solution.
-
Use the plotted time dependence to estimate the peak’s propagation velocity
$c$ . Compare the deduced$c$ to (21.4). -
Our solution of the wave equation for a plucked string leads to the formation of a wave packet that corresponds to the sum of multiple normal modes of the string. On the right in Figure 21.3 we show the motion resulting from the string initially placed in a single normal mode (standing wave),
$$\tag*{21.25} y(x,t=0) = 0.001 \sin 2\pi x, \quad\frac{\partial y}{\partial t}(x,t=0) =0,$$ and with friction (to be discussed soon) included.
Modify your program to incorporate this initial condition and see if a normal mode results.
-
Observe the motion of the wave for initial conditions corresponding to the sum of two adjacent normal modes. Does beating occur?
-
When a string is plucked near its end, a pulse reflects off the ends and bounces back and forth. Change the initial conditions of the model program to one corresponding to a string plucked exactly in its middle and see if a traveling or a standing wave results.
-
Figure 21.4 shows the wave packets that result as a function of time for initial conditions corresponding to the double pluck indicated on the left in the figure. Verify that initial conditions of the form
$$\tag*{21.26} \frac{y(x,t=0)}{0.005} =\begin{cases} 0, & 0.0\leq x \leq 0.1 , \ 10x -1, & 0.1 \leq x \leq 0.2 ,\ -10x +3, & 0.2 \leq x \leq 0.3 , \ 0, & 0.3 \leq x \leq 0.7 , \ 10x -7 , & 0.7\leq x\leq 0.8 , \ -10x +9, & 0.8\leq x \leq 0.9 ,\ 0, & 0.9 \leq x\leq 1.0 \end{cases}$$
lead to this type of a repeating pattern. In particular, observe whether the pulses move or just oscillate up and down.
The string problem we have investigated so far can be handled by either a numerical or an analytic technique. We now wish to extend the theory to include some more realistic physics. These extensions have only numerical solutions.
Plucked strings do not vibrate forever because there is friction in the
real world. Consider again the element of a string between
where
In Figure 21.3 we show the resulting motion of a string plucked in the middle when friction is included. Observe how the initial pluck breaks up into waves traveling to the right and to the left that are reflected and inverted by the fixed ends. Because those parts of the wave with the higher velocity experience greater friction, the peak tends to be smoothed out the most as time progresses.
Exercise: Generalize the algorithm used to solve the wave equation
to now include friction and check if the wave’s behavior seems physical
(damps in time). Start with
We have derived the propagation velocity for waves on a string as
To derive the equation for wave motion with variable density and
tension, consider again the element of a string (Figure 21.1 right) used
in our derivation of the wave equation. If we do not assume the tension
If
In §21.4.1 we will solve for the tension in a string as a result of gravity. Readers interested in an alternate easier problem that still shows the new physics may assume that the density and tension are proportional:
While we would expect the tension to be greater in regions of higher density (more mass to move and support), being proportional is clearly just an approximation. Substitution of these relations into (21.31) yields the new wave equation:
Here
Up until this point we have been ignoring the effect of gravity upon our
string’s shape and tension. This is a good approximation if there is
very little sag in the string, as might happen if the tension is very
high and the string is light. Even if there is some sag, our solution
for
Figure 21.5 Top: A uniform string suspended from its ends in a gravitational field assumes a catenary shape. Bottom: A force diagram of a section of the catenary at its lowest point. Because the ends of the string must support the entire weight of the string, the tension now varies along the string.
Consider a string of uniform density
The trick is to convert (21.35) to a differential equation that we can solve. We
do that by replacing the slope
Yet because
where
Here we have chosen the
It is this variation in tension that causes the wave velocity to change for different positions on the string.
Movie of wave motion of a plucked catenary with friction.
Figure 21.6 The wave motion of a plucked catenary with friction. (Courtesy of Juan Vanegas.)
We have given you the program EqStringAnimate.py
(Listing 21.1) that
solves the wave equation. Modify it to produce waves on a catenary
including friction for the assumed density and tension given by (21.32)
with CatFriction.py
and CatString.py
that do this.)
-
Look for some interesting cases and create surface plots of the results.
-
Describe in words how the waves dampen and how a wave’s velocity appears to change.
-
Normal modes: Search for normal-mode solutions of the variable-tension wave equation, that is, solutions that vary as
$$\tag*{21.42} u(x,t) = A \cos(\omega t)\sin(\gamma x).$$ Try using this form to start your program and see if you can find standing waves. Use large values for
$\omega$ . -
When conducting physics demonstrations, we set up standing-wave patterns by driving one end of the string periodically. Try doing the same with your program; that is, build into your code the condition that for all times
$$\tag*{21.43} y(x=0,t) = A \sin\omega t.$$ Try to vary
$A$ and$\omega$ until a normal mode (standing wave) is obtained. -
(For the exponential density case.) If you were able to find standing waves, then verify that this string acts like a high-frequency filter, that is, that there is a frequency below which no waves occur.
-
For the catenary problem, plot your results showing both the disturbance
$u(x,t)$ about the catenary and the actual height$y(x,t)$ above the horizontal for a plucked string initial condition. -
Try the first two normal modes for a uniform string as the initial conditions for the catenary. These should be close to, but not exactly, normal modes.
-
We derived the normal modes for a uniform string after assuming that $ k(x) = {\omega} /{c(x)}$ is a constant. For a catenary without too much
$x$ variation in the tension, we should be able to make the approximation$$\tag*{21.44} c(x)^2 \simeq \frac{T(x)} {\rho} = \frac{T_0\cosh ({x}/{d})}{\rho}.$$ See if you get a better representation of the first two normal modes if you include some
$x$ dependence in$k$ .
Problem: An elastic membrane is stretched across the top of a square
box of sides
where
The description of wave motion on a membrane is basically the same as
that of 1-D waves on a string discussed in § 21.2, only now we have wave
propagation in two directions. Consider Figure 21.7 showing a square
section of the membrane under tension
Figure 21.7 A small part of an oscillating membrane and the forces that act on it.
Although the tension is constant over the small area in Figure 21.7,
there will be a net vertical force on the segment if the angle of
incline of the membrane varies as we move through space. Accordingly,
the net force on the membrane in the
where
$$\begin{align} \tag*{21.47} \sin \theta \approx\ \tan \theta &= \frac{\partial u}{ \partial y}\Bigr |{y+\Delta y},\quad \sin \phi \approx\ \tan \phi = \frac{\partial u}{ \partial y}\Bigr |{y} , \ \Rightarrow \sum F_z(x_{fixed}) & = T\Delta x \Bigl ( \frac{\partial u}{ \partial y}\Bigr |{y+\Delta y}- \frac{\partial u}{ \partial y}\Bigr |{y} \Bigr ) \ \approx\ T\Delta x \frac{\partial^2 u}{\partial y^2} \Delta y.\tag*{21.48} \end{align}$$
Similarly, the net force in the
$$\tag*{21.49} \sum F_z (y_{fixed}) \ =\ T\Delta y \Bigl ( \frac{\partial u}{ \partial x}\Bigr |{x+\Delta x}- \frac{\partial u}{ \partial x}\Bigr |{x} \Bigr ) \ \approx\ T\Delta y \frac{\partial^2 u}{\partial x^2} \Delta x.$$
The membrane section has mass
This is the 2-D version of the wave equation (21.4) that we studied
previously in one dimension. Here
The analytic or numerical solution of the partial differential equation
(21.51) requires us to know both the boundary conditions and the initial
conditions. The boundary conditions hold for all times and were given
when we were told that the membrane is attached securely to a square box
of side
As required for a second order equation, the initial conditions has two
parts, the shape of the membrane at time
Second, we are told that the membrane is released from rest, which means:
where we write partial derivative because there are also spatial variations.
The analytic solution is based on the guess that because the wave
equation (21.51) has separate derivatives with respect to each
coordinate and time, the full solution
After substituting this into (21.51) and dividing by
The only way that the LHS of (21.57) can be true for all time while the RHS is also true for all coordinates, is if both sides are constant:
In (21.59) and (21.60) we have included the further deduction that
because each term on the RHS of (21.58) depends on either
We now apply the boundary conditions:
The fixed values for the eigenvalues
\begin{equation}\tag*{21.65} \xi^2 = q^2 + k^2 \hspace{6ex} \Rightarrow \hspace{6ex} \xi_{kq}=\pi\sqrt{k^2 +q^2}. \end{equation}
The full space-time solution now takes the form
\begin{equation}\tag*{21.66} u_{kq}= \left[G_{kq}\cos c\xi t +H_{kq} \sin c\xi t\right] \sin k x \sin q y, \end{equation}
where
\begin{equation}\tag*{21.67} u(x,y,t)=\sum_{k=1}^\infty \sum_{q=1}^\infty \left[G_{kq}\cos c\xi t +H_{kq} \sin c\xi t \right] \sin k x \sin q y. \end{equation}
While an infinite series is not a good algorithm, the
initial and boundary conditions means that only the
\begin{equation}\tag*{21.68} u(x,y,t)= \cos c\sqrt{5} \ \sin 2x\ \sin y, \end{equation}
where
**Figure 21.8 ** The standing wave pattern on a square box top at three different times.
Listing 21.1 Waves2D.py solves the wave equation numerically for a vibrating membrane.
The development of an algorithm for the solution of the 2-D wave equation (21.51) follows that of the 1-D equation in \S~21.2.2. We start by expressing the second derivatives in terms of central differences:
After discretizing the variables,
$$\tag*{21.72} \boxed{u^{k+1}{i,j} = 2u^k{i,j}-u^{k-1}{i,j} \frac{c^2}{c'^2} \Bigl [u^{k}{i+1,j}
- u^{k}{i-1,j} -4u^k{i,j}
- {u^{k}{i,j+1} + u^{k}{i,j-1} } \Bigr ],}$$
where as before
$$\tag*{21.73} 0 = \frac{\partial u (t=0)}{\partial t} \approx \frac{u^1_{i,j}-u^{-1}{i,j}} { 2\Delta t}, \hspace{6ex} \Rightarrow\hspace{6ex} u^{-1}{i,j} =u^1_{i,j}.$$
After substitution into (21.72) and solving for
$$\tag*{21.74} u^1_{i,j} = u^0_{i,j}+ \frac{c^2}{ 2 c'^2} \Bigl [u^{0}_{i+1,j}
- u^{0}{i-1,j} -4u^0{i,j}
- {u^{0}{i,j+1} + u^{0}{i,j-1} } \Bigr ].$$
Because the displacement
The program Wave2D.py
in Listing 21.2 solves the 2-D wave equation
using the time-stepping (leapfrog) algorithm. The program
Waves2Danal.py
computes the analytic solution. The shape of the
membrane at three different times are shown in Figure 21.8.
### VibratingMembrane.py, Notebook Version
import numpy as np
import matplotlib.pyplot as plt
import mpl_toolkits.mplot3d.axes3d
t=0 # initial time
den=390.0 #density kg/m2
ten=180.0 # tension N/m2
c=np.sqrt(ten/den) # propagation velocity
s5=np.sqrt(5)
N=32
t=0
def membrane(t,X,Y): # analytic solution vibrating membrane
return np.cos(c * s5* t) * np.sin( 2* X)*np.sin(Y)
plt.ion() # interactive on
fig=plt.figure()
ax = fig.add_subplot(111, projection='3d')
xs=np.linspace(0, np.pi,32) # 0 to pi increm: pi/32
ys=np.linspace(0, np.pi,32) # same
X, Y = np.meshgrid(xs,ys) # define x,grid
Z = membrane(0, X, Y) # initial condition t=0
wframe = None
ax.set_xlabel('X') # name x axis
ax.set_ylabel('Y')
ax.set_title('Vibrating membrane')
for t in np.linspace(0,50,200): #total time 10 divided by 40
oldcol=wframe # previous frame
Z=membrane(t,X,Y) # membrane at t !=0
wframe=ax.plot_wireframe(X,Y,Z)# plot the wireframe
# Remove old frame before drawing
if oldcol is not None:
ax.collections.remove(oldcol)
plt.draw() #plot new frame
print("finished")