|
Research 
|
Profile 
|
Teaching 
|
Software 
|
Contact 
|

[]

INSTITUTE OF MATHEMATICS


THE WAVE EQUATION

MY STUDENTS DATABASE

ashay@dharwadker.org

Copyright © by Ashay Dharwadker. All rights reserved.
 


 

The Wave Equation

Today's Link :: Friedrich Wilhelm Bessel

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Δy2u/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 = c22u.
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+...
v = 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.