Vibrating Drumhead: Symmetric Case

This notebook presents a solution to the wave equation which models a vibrating round drumhead that relies on Mathematica to separate, solve differential equations, and plot functions. We should be able to do the vast majority of this work by hand, if we like.

Clear[u,T,R]
Off[General::"spell1"]

Here is the wave equation in polar coordinates:

waveeq[u_]:=(1/(c^2))D[u,{t,2}]==D[u,{r,2}]+(1/r)D[u,{r,1}]+(1/(r^2))D[u,{θ,2}]

waveeq[u[r,θ,t]]

Solution: Symmetric Case

Separating Variables

If the drumhead is spherically symmetric, this simplifies since u(r,θ,t) = u(r,t). Let's separate variables,  u(r,t)=R(r)T(t)

waveeq[R[r]T[t]]

To separate, we always divide by u, so let's define an equation that does that for us.

waveeqsep[u_]:=(1/(c^2))((D[u,{t,2}])/u)==((D[u,{r,2}]+(1/r)D[u,{r,1}]+(1/(r^2))D[u,{θ,2}])/u)

sep1=Simplify[waveeqsep[R[r]T[t]]]

Now we have something we can work with. Since the left hand side of the equation depends only on t, and the right hand side depends only on r, they must be a constant.Let's choose that constant to be k.

deq1=Part[sep1,1]==k
deq2=Part[sep1,2]==k

We can use Mathematica to help us solve these differential equations. Let's start with the first.

sol1=DSolve[deq1,T[t],t]

Let's think about what this solution is telling us. If the constant k>0, then we will get exponential decay or exponential growth. Physically, we have a drumhead that is going to oscillate back and forth. To get that oscillation, we require k<0, which we can ensure if we set k=-(λ^2). Let's go back and solve the differential equations again, with this new constant (we go right back to the beginning to ensure Mathematica gets sines and cosines in its solution):

deq1b=Part[sep1,1]==-(λ^2)
deq2b=Part[sep1,2]==-(λ^2)

sol1=DSolve[deq1b,T[t],t]
T[t_]=Part[sol1,1,1,2]

OK, now we can tackle the second differential equation. It is much more complicated to solve by hand. Notice I use different constants in this solution, so we will not simply get C[1] and C[2] again.

sol2=DSolve[deq2b,R[r],r,GeneratedParameters->W]
R[r_]=Part[sol2,1,1,2]

There are those Bessel functions! We should be very happy and thankful that Mathematica has saved us the trouble of solving this differential equation, which would require (from first principles) a series solution which would lead to Bessel functions.

Determining the Constants of Integration

Now, let's think about this solution. The solution should remain bounded for r~0. Let's look at those Bessel functions more closely.

Plot[{BesselJ[0,r λ]/.λ->1,BesselY[0,r λ]/.λ->1},{r,0,1},PlotStyle->{Red,Blue},AxesLabel->{"r","J (red) and Y (blue)"},ImageSize->{500,Automatic}]

We see that the Bessel function (Subscript[Y, 0]) is unbounded at r=0; this leads us to eliminate it by choosing W[2]=0.

R[r_]=R[r]/.W[2]->0

Since we are looking for u(r,t)=R(r)T(t), that constant W[1] can be combined with other constants. We can therefore set W[1]=1 without loss of generality.

R[r_]=R[r]/.W[1]->1

Now, we need to consider what our boundary conditions are that will determine the constant of integration in the problem, and if there is anything else we should require (like earlier when we said k=-(λ^2)). Let's choose

u(a,t)=0, t≥0 (clamped circumference)
u(r,0)=f(r) (initial conditions are radially symmetric, depend only on r)
(∂/(∂t))u(r,0)=g(r)

Let's consider that first equation.

u(a,t)=0
T(t) R(a) = 0
R(a)=0 (if the above is to be true for all values of t)

This has an impact on the R(r) we found above. Let's sketch the Bessel function again:

Plot[BesselJ[0,x],{x,0,15},PlotStyle->{Red,Blue},AxesLabel->{"x","(Subscript[J, 0])(x)"},ImageSize->{500,Automatic}]

We see that (Subscript[J, 0])(x)has infinitely many zeros. Therefore, the boundary of the drum must occur at one of the zeros to satisfy the condition u(a,t)=0.The roots of the Bessel function are:

Solve[BesselJ[0,x]==0,x]

hard to find. Let's call the roots of the Bessel function (Subscript[α, i]), where (Subscript[α, i])< (Subscript[α, i+1]), i=1,2,3,4,....

Therefore, if

a λ=(Subscript[α, n]) -> (Subscript[λ, n])=(Subscript[α, n])/a, where (Subscript[α, n]) are the zeros of the Bessel function (Subscript[J, 0]).

then the boundary condition that the circumference is clamped will be satisfied.This leads to a "quantization", or multiple solutions to the problem, for different values of the constant n:

(Subscript[R, n])[r_]=R[r]/.λ->(Subscript[λ, n])

OK, so all that is left is to determine the constant in T(t). With the new requirements on λ, T(t) becomes:

(Subscript[T, n])[t_]=T[t]/.λ->(Subscript[λ, n])

and with each different value of n, we have a different function. Therefore, the constants should depend on n as well:

(Subscript[T, n])[t_]=(Subscript[T, n])[t]/.C[1]->(Subscript[A, n]) /.C[2]->(Subscript[B, n])

Let's combine to get our (Subscript[u, n])(r,t), which is

(Subscript[u, n])[r_,t_]=(Subscript[R, n])[r](Subscript[T, n])[t]

The complete solution is the sum over all the different values of n which are possible:

u[r_,t_]=(∑_(n=1)^∞)(Subscript[u, n])[r,t]

Now, we need to consider what our boundary conditions will do to the solution. Let's tackle them in turn.

u(r,0)=f(r)

beq1=u[r,0]==f[r]

This is similar to a Fourier sine series, except it involves Bessel functions instead of sines. This analogy gives us the technique we can use to determine (Subscript[A, n]).

((∑_(n=1)^∞)(Subscript[J, 0])(r (Subscript[λ, n])) (Subscript[A, n])==f[r])*(Subscript[J, 0])(r (Subscript[λ, m])) 
(∑_(n=1)^∞)(Subscript[J, 0])(r (Subscript[λ, n])) (Subscript[A, n])(Subscript[J, 0])(r (Subscript[λ, m])) ==f[r]*(Subscript[J, 0])(r (Subscript[λ, m])) 
(∫_0^a)(∑_(n=1)^∞)(Subscript[A, n])*(Subscript[J, 0])(r (Subscript[λ, n]))(Subscript[J, 0])(r (Subscript[λ, m]))dr ==(∫_0^a)f[r]*(Subscript[J, 0])(r (Subscript[λ, m])) dr
(∑_(n=1)^∞)(Subscript[A, n])*(∫_0^a)(Subscript[J, 0])(r (Subscript[λ, n])) (Subscript[J, 0])(r (Subscript[λ, m]))dr ==(∫_0^a)f[r]*(Subscript[J, 0])(r (Subscript[λ, m])) dr

(∫_0^a)BesselJ[0,r (Subscript[λ, n])] BesselJ[0,r (Subscript[λ, m])]\[DifferentialD]r

Looks like Mathematica can't do this integral. Let's back up and try again. A quick look at the properties of Bessel functions (see http://mathworld.wolfram.com/BesselFunctionoftheFirstKind.html Eq. (51)) leads us to realize that we should include an r, and then exploit the orthogonality properties of the Bessel function:

((∑_(n=1)^∞)(Subscript[J, 0])(r (Subscript[λ, n])) (Subscript[A, n])==f(r))*(Subscript[J, 0])(r (Subscript[λ, m])) r
(∑_(n=1)^∞)(Subscript[J, 0])(r (Subscript[λ, n])) (Subscript[A, n])(Subscript[J, 0])(r (Subscript[λ, m])) r==f(r)*(Subscript[J, 0])(r (Subscript[λ, m])) r
(∫_0^a)(∑_(n=1)^∞)(Subscript[A, n])*(Subscript[J, 0])(r (Subscript[λ, n])) (Subscript[J, 0])(r (Subscript[λ, m])) r dr ==(∫_0^a)f(r)*(Subscript[J, 0])(r (Subscript[λ, m]))r dr
(∑_(n=1)^∞)(Subscript[A, n])*(∫_0^a)(Subscript[J, 0])(r (Subscript[λ, n])) (Subscript[J, 0])(r (Subscript[λ, m])) r dr ==(∫_0^a)f(r)*(Subscript[J, 0])(r (Subscript[λ, m])) rdr
(∑_(n=1)^∞)(Subscript[A, n])* ((a^2)/2)((((Subscript[J, 0])(r (Subscript[λ, n]))))^2) (Subscript[δ, m,n]) ==(∫_0^a)f(r)*(Subscript[J, 0])(r (Subscript[λ, m]))rdr
(Subscript[A, m])* ((a^2)/2)((((Subscript[J, 1])(r (Subscript[λ, m]))))^2) ==(∫_0^a)f(r)*(Subscript[J, 0])(r (Subscript[λ, m]))rdr
(Subscript[A, m])==(2/( (a^2)((((Subscript[J, 1])(r (Subscript[λ, m]))))^2) ))(∫_0^a)f(r)*(Subscript[J, 0])(r (Subscript[λ, m]))rdr

You can do exactly the same procedure for the other boundary condition:

(d/(d t))u(r,0)=g(r)

beq2=D[u[r,t],{t,1}]==g[r]/.t->0

Proceeding as above, we determine

(Subscript[B, m])==(2/( c  (Subscript[λ, m])(a^2)((((Subscript[J, 1])(r (Subscript[λ, m]))))^2) ))(∫_0^a)g(r)*(Subscript[J, 0])(r (Subscript[λ, m]))rdr

We have solved the problem!

u[r,t]
(Subscript[A, m_]):=(2/( (a^2)((BesselJ[1,a (Subscript[λ, m])])^2) ))(∫_0^a)f[r]*BesselJ[0,r (Subscript[λ, m])]r\[DifferentialD]r
(Subscript[B, m_]):=(2/( c  (Subscript[λ, m])(a^2)((BesselJ[1,a (Subscript[λ, m])])^2) ))(∫_0^a)g[r]*BesselJ[0,r (Subscript[λ, m])]r\[DifferentialD]r

Actually doing these integrals can be complicated.

Example 1

Assume a circular drum head, clamped at the edges, is initially flat, with initial velocity of -100m/s due to an explosion which takes place near the drumhead. Take a=1, c=100, and determine how the membrane vibrates.

We need the zeros of the Bessel function, which we saw earlier gave Mathematica a headache. But they are so important that Mathematica has a way of accessing them.

a=1;
c=100;
f[r_]=0;
g[r_]=-100;
(Subscript[λ, m_]):=((BesselJZero[0,m])/a)//N
(Subscript[A, m_]):=(2/( (a^2)((BesselJ[1,a (Subscript[λ, m])])^2) ))(∫_0^a)f[r]*BesselJ[0,r (Subscript[λ, m])]r\[DifferentialD]r
(Subscript[B, m_]):=(2/( c  (Subscript[λ, m])(a^2)((BesselJ[1,a (Subscript[λ, m])])^2) ))(∫_0^a)g[r]*BesselJ[0,r (Subscript[λ, m])]r\[DifferentialD]r
u[r,t]:=(∑_(n=1)^∞)BesselJ[0,r (Subscript[λ, n])] (Cos[c t (Subscript[λ, n])] (Subscript[A, n])+Sin[c t (Subscript[λ, n])] (Subscript[B, n]))

Let's have Mathematica do the integrals for us. We have to truncate the series.

sol[r_,t_]=(∑_(n=1)^10)BesselJ[0,r (Subscript[λ, n])] (Cos[c t (Subscript[λ, n])] (Subscript[A, n])+Sin[c t (Subscript[λ, n])] (Subscript[B, n]))

Is this oscillation periodic? Let's explore:

Plot[sol[0.2,t],{t,0,0.2}]

There is an oscillation, but it does not appear to be periodic. The following arithmetic seems to bear this result out:

2Pi/240.4825557695773
4Pi/552.007811028631
8Pi/865.3727912911013

Plots of the vibrating drumhead:

Do[

plot1[t]=RevolutionPlot3D[sol[r,t],{r,0,1},{phi,0,2Pi},PlotRange->{{-1,1},{-1,1},{-1.5,1.5}},Boxed->False,Axes->False,ImageSize->{400,Automatic}],

{t,0,5,5/20}]

Manipulate[Show[plot1[t]],{t,0,5,5/20}]

Example 2

Assume a circular drum head, clamped at the edges, which is initially flat. The initial velocity of the drumhead is given by

g(r) ={({{1, if 0<r<(1/2)}, {0, if (1/2)<r<1}})

Take a=1, c=1, and determine how the membrane vibrates.

a=1;
c=1;
f[r_]=0;
g[r_]=If[0<r<(1/2),1,0];
(Subscript[λ, m_]):=((BesselJZero[0,m])/a)//N
(Subscript[A, m_]):=(2/( (a^2)((BesselJ[1,a (Subscript[λ, m])])^2) ))(∫_0^a)f[r]*BesselJ[0,r (Subscript[λ, m])]r\[DifferentialD]r
(Subscript[B, m_]):=(2/( c  (Subscript[λ, m])(a^2)((BesselJ[1,a (Subscript[λ, m])])^2) ))(∫_0^a)g[r]*BesselJ[0,r (Subscript[λ, m])]r\[DifferentialD]r
u[r,t]:=(∑_(n=1)^∞)BesselJ[0,r (Subscript[λ, n])] (Cos[c t (Subscript[λ, n])] (Subscript[A, n])+Sin[c t (Subscript[λ, n])] (Subscript[B, n]))

Let's have Mathematica do the integrals for us. We have to truncate the series.

sol[r_,t_]=(∑_(n=1)^10)BesselJ[0,r (Subscript[λ, n])] (Cos[c t (Subscript[λ, n])] (Subscript[A, n])+Sin[c t (Subscript[λ, n])] (Subscript[B, n]))

Is this oscillation periodic? Let's explore:

Plot[sol[0,t],{t,0,20}]

There is an oscillation, but it does not appear to be periodic.

Do[

plot1[t]=RevolutionPlot3D[sol[r,t],{r,0,1},{phi,0,2Pi},PlotRange->{{-1,1},{-1,1},{-0.5,0.5}},Boxed->False,Axes->False,ImageSize->{400,Automatic}],

{t,0,5,5/20}]

Manipulate[Show[plot1[t]],{t,0,5,5/20}]

Example 3

Assume a circular drum head, clamped at the edges, which is initially displaced by the function f(r)=1-(r^2). The initial velocity of the drumhead is zero.

Take a=1, c=1, and determine how the membrane vibrates.

a=1;
c=1;
f[r_]=1-(r^2);
g[r_]=0;
(Subscript[λ, m_]):=((BesselJZero[0,m])/a)//N
(Subscript[A, m_]):=(2/( (a^2)((BesselJ[1,a (Subscript[λ, m])])^2) ))(∫_0^a)f[r]*BesselJ[0,r (Subscript[λ, m])]r\[DifferentialD]r
(Subscript[B, m_]):=(2/( c  (Subscript[λ, m])(a^2)((BesselJ[1,a (Subscript[λ, m])])^2) ))(∫_0^a)g[r]*BesselJ[0,r (Subscript[λ, m])]r\[DifferentialD]r
u[r,t]:=(∑_(n=1)^∞)BesselJ[0,r (Subscript[λ, n])] (Cos[c t (Subscript[λ, n])] (Subscript[A, n])+Sin[c t (Subscript[λ, n])] (Subscript[B, n]))

Let's have Mathematica do the integrals for us. We have to truncate the series.

sol[r_,t_]=(∑_(n=1)^10)BesselJ[0,r (Subscript[λ, n])] (Cos[c t (Subscript[λ, n])] (Subscript[A, n])+Sin[c t (Subscript[λ, n])] (Subscript[B, n]))

Is this oscillation periodic? Let's explore:

Plot[sol[0,t],{t,0,10}]

There is an oscillation, but it does not appear to be periodic.

Do[

plot1[t]=RevolutionPlot3D[sol[r,t],{r,0,1},{phi,0,2Pi},PlotRange->{{-1,1},{-1,1},{-1.5,1.5}},Boxed->False,Axes->False,ImageSize->{400,Automatic}],

{t,0,5,5/20}]

Manipulate[Show[plot1[t]],{t,0,5,5/20}]