The Wave Equation
The Problem
We study the motion of a
stretched circular membrane, such as the vibrating surface of a drum, subject
to the following physical assumptions:
-
The membrane has uniform mass
per unit area, is perfectly flexible and so thin that it does not offer
any resistance to bending.
-
The membrane is stretched and
then fixed along its circular boundary in the xy-plane, the tension T per
unit length being the same at all points and along all directions, not
changing during the motion.
-
The vertical deflection u = u(x,y,t)
of the point (x,y) of the membrane at time t is small compared to the diameter.
The aim is to find an explicit
formula for u(x,y,t), given the initial deflection, so that we may model
the motion of the vibrating drum head.
|
The Differential
Equation
| Consider a small portion
of the membrane:
The sides of the portion
are of small lengths Δx
and Δy.
The tension T is the force per unit length, hence the forces acting on
the edges are TΔx
and TΔy,
tangent to the membrane since it is perfectly flexible. First consider
the components of the forces parallel to the xy-plane. Along the x-axis,
the resultant TΔy
cos(β) -
TΔy cos(α)
is (almost) zero, since α
and β are
small. Similarly, the resultant force is (almost) zero along the y-axis.
Thus, the motion of the portion is completely governed by the vertical
components of the forces parallel to the u-axis. The vertical resultant
force for the edges parallel to the x-axis is
TΔy
sin(β) -
TΔy sin(α)
= TΔy(sin(β)
- sin(α))
= TΔy(tan(β)
- tan(α))
(since α
and β are
small)
= TΔy(ux(x +Δx,y1) -ux(x,y2))
(by Rolle's Mean Value Theorem)
where ux denotes
the partial derivative of u with respect to x and y1 and y2
are values between y and y+Δy.
Similarly, the vertical resultant force for the edges parallel to the y-axis
is
TΔx(uy(x1,y+Δy)-uy(x2,y))
where uy denotes
the partial derivative of u with respect to y and x1 and x2
are values between x and x+Δx.
Thus the total resultant force on the portion is given by
F = TΔy(ux(x+Δx,y1)-ux(x,y2))+TΔx(uy(x1,y+Δy)-uy(x2,y)).
Let ρ
denote the mass per unit area of the membrane. The mass of the portion
is m = ρΔxΔy
and its acceleration is a = ∂2u/∂t2.
By Newton's second law of motion F = ma we have
ρΔxΔy∂2u/∂t2 = TΔy(ux(x+Δx,y1)-ux(x,y2))+TΔx(uy(x1,y+Δy)-uy(x2,y)).
Dividing both sides by ρΔxΔy
we get
∂2u/∂t2 = (T/ρ)((ux(x+Δx,y1)-ux(x,y2))/Δx)+(uy(x1,y+Δy)-uy(x2,y))/Δy).
Now let Δx
and Δy approach
zero to obtain the differential equation
∂2u/∂t2 = c2(∂2u/∂x2+∂2u/∂y2)
where c2 = T/ρ.
This differential equation is
called the 2-dimensional wave equation. |
The Laplacian
The operator ∇ = ∂2/∂x2+∂2/∂y2
is called the Laplacian and the 2-dimensional wave equation may
then be written in terms of the Laplacian simply as
∂2u/∂t2 = c2∇2u.
Since the boundary of the circular
membrane may be expressed by the simple equation r = constant, we first transform
the Laplacian into polar coordinates (r,θ).
Note that x = r cos(θ)
and y = r sin(θ)
and we denote partial derivatives by subscripts.
By the chain rule,
ux = urrx+uθθx.
Differentiating again using
the product rule,
(1)
uxx = urxrx+urrxx+uθxθx+uθθxx.
Again by the chain rule,
urx = urrrx +urθθx
and uθx = uθrrx+uθθθx.
Now r = (x2+y2)1/2
and θ = tan-1(y/x),
so rx = x/r and θx = -y/r2.
Differentiating again, rxx = (r2-x2)/r3 = y2/r3
and θxx = 2xy/r4.
Substituting into (1) and assuming the continuity of first and second derivatives
so that urθ = uθr,
we obtain
(2)
uxx = (x2/r2)urr-2(xy/r3)urθ+(y2/r4)uθθ+(y2/r3)ur +2(xy/r4)uθ.
Similarly,
(3)
uyy = (y2/r2)urr+2(xy/r3)urθ+(x2/r4)uθθ+(x2/r3)ur -2(xy/r4)uθ.
Now add (2) and (3) to find
the Laplacian of u in polar coordinates:
∇2u = uxx+uyy = urr+(1/r)ur+(1/r2)uθθ
∇2u = ∂2u/∂r2+(1/r)∂u/∂r+(1/r2)∂2u/∂θ2
|
The Solution
Using polar coordinates
(r,θ) where
x = r cos(θ)
and y = r sin(θ),
the wave equation becomes
∂2u/∂t2 = c2(∂2u/∂r2+(1/r)∂u/∂r+(1/r2)∂2u/∂θ2).
Consider solutions of the wave
equation that are radially symmetric, i.e. do not depend on θ.
Then the wave equation becomes
(4)
∂2u/∂t2 = c2(∂2u/∂r2+(1/r)∂u/∂r)
and we wish to find solutions
u(r,t). Since the membrane is fixed along the boundary r = 1, we have the
boundary condition
(5)
u(1,t) = 0 for all t ≥ 0.
Solutions that do not depend
on θ will
occur if the initial conditions do not depend on θ
(6)
u(r,0) = f(r) (Initial Deflection)
(7)
ut(0) = g(r) (Initial Velocity)
Using the method of separating
variables, we first find solutions of (4) that satisfy the boundary condition
(5). Write
u(r,t) = v(r)w(t)
Now ut = vwt,
utt = vwtt, ur = vrw, urr = vrrw
so that (4) becomes
vwtt = c2(vrrw+(1/r)vrw)
vwtt = c2w(vrr+(1/r)vr)
(wtt/(c2w)) = (1/v)(vrr+(1/r)vr)
separating the variables. The
expressions on both sides must be equal to a negative constant -k2
in order to obtain solutions that satisfy the boundary condition (5) without
being identically zero. Hence
(wtt/(c2w)) = (1/v)(vrr+(1/r)vr) = -k2.
We now have two ordinary differential
equations
(8)
wtt+λ2w = 0
where λ = ck
(9)
vrr+(1/r)vr+k2v = 0
First solve (9) as follows.
Define a new variable s = kr. Then by the chain rule vr = vssr = vsk
and vrr = vssk2. Substituting in (9)
vssk2+(1/(s/k))vsk+k2v = 0
(10) vss+(1/s)vs+v = 0
we obtain Bessel's Differential
Equation (10). This equation can be solved using the Frobenius power
series method. Assume that there is a solution of the form
v
= c0+c1s1+c2s2+c3s3+c4s4+...
vs = c1+2c2s1+3c3s2+4c4s3+...
vss = 2c2+2.3c3s1+3.4c4s2+...
Substituting in (10)
(2c2+2.3c3s1+3.4c4s2+...)+(1/s)(c1+2c2s1+3c3s2+...)+(c0+c1s1+c2s2+...) = 0
c1s-1+(2.2c2+c0)s0+(3.3c3+c1)s1+(4.4c4+c2)s2+... = 0
Equating the coefficients to
zero
c1 = 0
2.2c2+c0 = 0
3.3c3+c1 = 0
4.4c4+c2 = 0
.
.
.
Now c0 is arbitrary,
so put c0 = 1. Hence
c0 = 1
c1 = 0
c2 = -1/(2.2) = -1/(221!2)
c3 = 0
c4 = 1/(2.2.4.4) = 1/(242!2)
c5 = 0
c6 = -1/(2.2.4.4.6.6) = -1/(263!2)
.
.
.
In general, for n = 0,1,2,3,...
c2n = (-1)n/(22nn!2)
c2n+1 = 0
Thus, we have found a solution
of Bessel's equation (10) known as a Bessel Function of the First Kind
and traditionally denoted by
J0(s) = v(s) = n = 0Σ∞((-1)n/(22nn!2))s2n.
Hence we have found a solution
to the differential equation (9), substituting s = kr, given by
(11)
v(r) = J0(kr).
Now, the boundary condition
(5) implies
u(1,t) = 0 for all
t ≥ 0
v(1)w(t) = 0 for all t ≥ 0.
Since w(t) = 0 for all t ≥ 0
would imply u(r,t) = 0 or all t ≥ 0,
we require v(1) = 0. Hence v(1) = J0(k) = 0. There are infinitely
many roots of this equation, for example, k = 2.4048, 5.5201, 8.6537, 11.7915,
14.9309,... correct to four decimal places.
We now solve (8) for any
of the above values of k, as follows. Suppose w(t) = eqt is a
solution of (8) for some suitable choice of q. We have wt = qeqt
and wtt = q2eqt. Substituting in (8) we
have
q2eqt+c2k2eqt = 0
(q2+c2k2)eqt = 0.
Hence q2+c2k2 = 0.
This equation has two complex conjugate roots q = cki and q = -cki.
Thus, using Euler's formula, we have two complex solutions
w1(t) = ei
ckt = cos(ckt)+i
sin(ckt)
w2(t) = e-i
ckt = cos(ckt)-i
sin(ckt)
Now, using the linearity of
differentiation, w(t) = (1/2)(w1+w2)+(1/2i
)(w1-w2) = cos(ckt)+sin(ckt)
is a real solution of (8). In fact,
(12)
w(t) = a cos(ckt) + b sin(ckt)
is a solution of (8) for arbitrary
constants a and b. We have thus found a solution for the wave equation
(4) satisfying the boundary condition (5) from (11) and (12) by
u(r,t) = v(r)w(t) = J0(kr).(a
cos(ckt) + b sin(ckt)).
This solution u(r,t) is called
an eigenfunction of our problem and the corresponding λ = ck
is called an eigenvalue. The vibration of the membrane corresponding
to u(r,t) has the frequency λ/2π
cycles per unit time. |
The Model
The following C++ program
models the solution u(r,t) found above. We compute the Bessel functions
of the first kind and determine the constants a and b by solving the Bessel-Fourier
series numerically. We assume that the initial deflection (6) of the membrane
is
u(r,0) = f(r) = (1-r2)1/2
and the initial velocity (7)
of the membrane is
ut(0) = g(r) = 0.
The graphics are implemented
using OpenGL materials and shading with smooth normals. |
Bessel.cpp
#include <iostream>
#include <fstream>
#include <string>
#include <cmath>
using namespace std;
float factorial(float n);
float J(float n, float x);
float f(float r);
float g(float r);
float a(float m, float R);
float b(float m, float R);
float u(float m, float R,
float
r,
float
t);
ofstream outfile("data.txt");
int main()
{
for(float t = 0;
t<2.5; t+ = 0.1)
{
for(float r = 0;
r< = 3; r+ = 0.1)
{
outfile<<u(1,1,r,t)<<"
";
}
outfile<<endl;
}
system("PAUSE");
return
0;
}
float factorial(float
n)
{
float product = 1;
for(int i = 1;
i< = n; i++)
product* = i;
return product;
}
float J(float n, float x)
{
float sum = 0;
for(float m = 0;
m<10; m+ = 1)
sum+ = (pow(-1,m)*pow(x,2*m))/(pow(2,2*m+n)*factorial(m)*factorial(n+m));
return pow(x,n)*sum;
}
float f(float r)
{
return 1-pow(r,2);
}
float g(float r)
{
return 0;
}
float a(float m, float R)
{
float alpha;
if(m = 1)
alpha = 2.4048;
else if(m = 2)
alpha = 5.5201;
else if(m = 3)
alpha = 8.6537;
else if(m = 4)
alpha = 11.7915;
else if(m = 5)
alpha = 14.9309;
else alpha = 0;
float integral = 0;
float step = R/10;
for(float r = 0;
r< = R; r+ = step)
integral+ = r*f(r)*J(0,alpha*r/R)*step;
return (2/(pow(R,2)*pow(J(1,alpha),2)))*integral;
}
float b(float m, float R)
{
float c = 1;
float alpha;
if(m = 1)
alpha = 2.4048;
else if(m = 2)
alpha = 5.5201;
else if(m = 3)
alpha = 8.6537;
else if(m = 4)
alpha = 11.7915;
else if(m = 5)
alpha = 14.9309;
else alpha = 0;
float integral = 0;
float step = R/10;
for(float r = 0;
r< = R; r+ = step)
integral+ = r*g(r)*J(0,alpha*r/R)*step;
return (2/(c*alpha*R*pow(J(1,alpha),2)))*integral;
}
float u(float m, float R,
float
r,
float
t)
{
float c = 1;
float alpha;
if(m = 1)
alpha = 2.4048;
else if(m = 2)
alpha = 5.5201;
else if(m = 3)
alpha = 8.6537;
else if(m = 4)
alpha = 11.7915;
else if(m = 5)
alpha = 14.9309;
else alpha = 0;
float lambda = c*alpha/R;
return (a(m,R)*cos(lambda*t)+b(m,R)*sin(lambda*t))*J(0,lambda*r/c);
} |
Exercise
| A smaller drum should have a higher fundamental frequency than a larger drum,
tension and density being the same. How can we conclude this from our formulas?
Modify the program to demonstrate this visually for drums of various sizes. |
|
