Vibrating Drumhead: General 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: General Case

Separating Variables

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=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 and θ, they must be a constant. Let's choose that constant to be k=-(λ^2), since other values of the constant k will not produce the required periodic solutions (see discussion for symmetric case for details).

deq1=Part[sep1,1]==-(λ^2)
deq2=Part[sep1,2]==-(λ^2)

We can use Mathematica to help us solve this differential equation.

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

Our deq2 is separable again.

Simplify[deq2]

Notice that we can rewrite this as

(r^2) (λ^2)+((r R'[r]+(r^2) R''[r])/(R[r]))=-((Θ''[θ])/( Θ[θ]))==(Subscript[k, 2])

OK, let's see what happens when we solve the theta part:

DSolve[-((Θ''[θ])/( Θ[θ]))==(Subscript[k, 2]),Θ[θ],θ]

This is fine, and we see that we are getting an oscillatory solution, which is what we expect physically, as long as (Subscript[k, 2])>0. However, a small tweaking of the constant will lead to a more pleasing solution, and ensure positivity. Let (Subscript[k, 2])=(μ^2). We will revisit this below when we start looking at determining the constants of integrations that arise in the problem.

sol2 =DSolve[-((Θ''[θ])/( Θ[θ]))==(μ^2),Θ[θ],θ,GeneratedParameters->W]
Θ[θ_]=Part[sol2,1,1,2]

OK, now we can tackle the final ordinary differential equation. It is much more complicated to solve by hand than the first two.

sol3=DSolve[(r^2) (λ^2)+((r R'[r]+(r^2) R''[r])/(R[r]))==(μ^2),R[r],r,GeneratedParameters->Z]
R[r_]=Part[sol3,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. Notice that for the symmetric case we got Bessel functions of order 0; not surprisingly, here we get Bessel functions of order μ.

Here is the solution.

R[r]Θ[θ]T[t]

Determining the Constants of Integration

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 the following:

1. Bounded at r=0

Now, let's think about this solution. We have lots of constants that we need to determine based on the physical system or boundary conditions, or initial conditions! The solution should remain bounded for r~0. Let's look at those Bessel functions more closely.

Plot[{BesselJ[1,r λ]/.λ->1,BesselY[1,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, 1]) is unbounded at r=0; this is also true for (Subscript[Y, m])for any value of m; this leads us to eliminate it by choosing Z[2]=0.

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

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

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

2. Clamped Circumference

Let's consider the clamped circumference equation. It leads to the following:

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

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

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

We see that (Subscript[J, 2.3])(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[2.3,x]==0,x]

hard to find. Let's call the roots of the Bessel function (Subscript[J, μ])(x) (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, μ])(x),

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])
Θ[θ_]=Θ[θ]
(Subscript[T, μ_,n_])[t_]=T[t]/.λ->(Subscript[λ, μ,n])

Until now, I have left the μ as a real number.

3. Periodicity in θ

Next up is the peridocity in the θ.

u(r,0,t) = u(r,2π,t) 
R(r) Θ(0)T(t) =  R(r) Θ(2π)T(t)
 Θ(0) = Θ(2π)

Also, we have

(∂/(∂θ))u(r,0,t) = (∂/(∂θ))u(r,2π,t)
(∂/(∂θ))Θ(0) = (∂/(∂θ))Θ(2π)

This translates into the two equations:

eq1=Θ[0]==Θ[2π]
eq2=Derivative[1][Θ][0]==Derivative[1][Θ][2π]

These are satisfied if μ is an integer. Therefore, we need to let μ=m, m=0,1,2,3,4,.... Again, this "quantizes" the system in a sense. Now we have the following:

(Subscript[R, m,n])[r]
(Subscript[Θ, m_])[θ_]=Θ[θ]/.μ->m
(Subscript[T, m,n])[t]

4. Initial Conditions

Here is where we are so far:

(Subscript[u, m_,n_])[r_,θ_,t_]=(Subscript[R, m,n])[r](Subscript[Θ, m])[θ](Subscript[T, m,n])[t]

OK, so all that is left is to determine the constants C[1], C[2], W[1], and W[2], and with each different value of m,n, we have a different function. Therefore, the constants should depend on m,n as well. The complete solution is the sum over m,n:

(Subscript[u, m_,n_])[r_,θ_,t_]=(Subscript[u, m,n])[r,θ,t]/.C[1]->(Subscript[A, m,n])/.C[2]->(Subscript[B, m,n])/.W[1]->(Subscript[(OverBar[A]), m,n])/.W[2]->(Subscript[(OverBar[B]), m,n])

The complete solution is the sum over m,n:

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

Now, we need to consider what our boundary conditions will do to the solution.

u(r,θ,0)=f(r,θ)

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

The manipulations which follow are significantly complicated that I have them in a pdf file, for ease of viewing and understanding. The file can be downloaded here.

We have solved the problem!

u[r_,θ_,t_]:=(∑_(m=0)^∞)(∑_(n=1)^∞)BesselJ[m,r (Subscript[λ, m,n])]((Subscript[a, m,n])Cos[m θ]+(Subscript[b, m,n])Sin[m θ])Cos[c (Subscript[λ, m,n]) t]+(∑_(m=0)^∞)(∑_(n=1)^∞)BesselJ[m,r (Subscript[λ, m,n])]((Subscript[a2, m,n])Cos[m θ]+(Subscript[b2, m,n])Sin[m θ])Sin[c (Subscript[λ, m,n]) t]

(Subscript[a, 0,n_]):=(2/( π (a^2)((BesselJ[1,a (Subscript[λ, 0,n])])^2) ))(∫_0^a)(∫_0^(2π))f[r,θ]BesselJ[0,r (Subscript[λ, 0,n])]r\[DifferentialD]θ\[DifferentialD]r
(Subscript[a, m_,n_]):=(2/( π (a^2)((BesselJ[m+1,a (Subscript[λ, m,n])])^2) ))(∫_0^a)(∫_0^(2π))f[r,θ]Cos[m θ]BesselJ[m,r (Subscript[λ, m,n])]r\[DifferentialD]θ\[DifferentialD]r
(Subscript[b, m_,n_]):=(2/( π (a^2)((BesselJ[m+1,a (Subscript[λ, m,n])])^2) ))(∫_0^a)(∫_0^(2π))f[r,θ]Sin[m θ]BesselJ[m,r (Subscript[λ, m,n])]r\[DifferentialD]θ\[DifferentialD]r

(Subscript[a2, 0,n_]):=(2/(c (Subscript[λ, 0,n]) π (a^2)((BesselJ[1,a (Subscript[λ, 0,n])])^2) ))(∫_0^a)(∫_0^(2π))g[r,θ]BesselJ[0,r (Subscript[λ, 0,n])]r\[DifferentialD]θ\[DifferentialD]r
(Subscript[a2, m_,n_]):=(2/(c (Subscript[λ, m,n]) π (a^2)((BesselJ[m+1,a (Subscript[λ, m,n])])^2) ))(∫_0^a)(∫_0^(2π))g[r,θ]Cos[m θ]BesselJ[m,r (Subscript[λ, m,n])]r\[DifferentialD]θ\[DifferentialD]r
(Subscript[b2, m_,n_]):=(2/(c (Subscript[λ, m,n]) π (a^2)((BesselJ[m+1,a (Subscript[λ, m,n])])^2) ))(∫_0^a)(∫_0^(2π))g[r,θ]Sin[m θ]BesselJ[m,r (Subscript[λ, m,n])]r\[DifferentialD]θ\[DifferentialD]r

Actually doing these integrals can be complicated, but you have probably already guessed that!

Example 1

Solve the general vibrating drumhead problem when a=c=1,f(r,θ)=(1-(r^2))r sin(θ), g(r,θ)=(1-(r^2))(r^2)sin(2θ).

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=1;
f[r_,θ_]=(1-(r^2))r Sin[θ];
g[r_,θ_]=(1-(r^2))(r^2) Sin[2θ];

(Subscript[λ, m_,n_]):=((BesselJZero[m,n])/a)//N

(Subscript[a, 0,n_]):=(2/( π (a^2)((BesselJ[1,a (Subscript[λ, 0,n])])^2) ))(∫_0^a)(∫_0^(2π))f[r,θ]BesselJ[0,r (Subscript[λ, 0,n])]r\[DifferentialD]θ\[DifferentialD]r
(Subscript[a, m_,n_]):=(2/( π (a^2)((BesselJ[m+1,a (Subscript[λ, m,n])])^2) ))(∫_0^a)(∫_0^(2π))f[r,θ]Cos[m θ]BesselJ[m,r (Subscript[λ, m,n])]r\[DifferentialD]θ\[DifferentialD]r
(Subscript[b, m_,n_]):=(2/( π (a^2)((BesselJ[m+1,a (Subscript[λ, m,n])])^2) ))(∫_0^a)(∫_0^(2π))f[r,θ]Sin[m θ]BesselJ[m,r (Subscript[λ, m,n])]r\[DifferentialD]θ\[DifferentialD]r

(Subscript[a2, 0,n_]):=(2/(c (Subscript[λ, 0,n]) π (a^2)((BesselJ[1,a (Subscript[λ, 0,n])])^2) ))(∫_0^a)(∫_0^(2π))g[r,θ]BesselJ[0,r (Subscript[λ, 0,n])]r\[DifferentialD]θ\[DifferentialD]r
(Subscript[a2, m_,n_]):=(2/(c (Subscript[λ, m,n]) π (a^2)((BesselJ[m+1,a (Subscript[λ, m,n])])^2) ))(∫_0^a)(∫_0^(2π))g[r,θ]Cos[m θ]BesselJ[m,r (Subscript[λ, m,n])]r\[DifferentialD]θ\[DifferentialD]r
(Subscript[b2, m_,n_]):=(2/(c (Subscript[λ, m,n]) π (a^2)((BesselJ[m+1,a (Subscript[λ, m,n])])^2) ))(∫_0^a)(∫_0^(2π))g[r,θ]Sin[m θ]BesselJ[m,r (Subscript[λ, m,n])]r\[DifferentialD]θ\[DifferentialD]r

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

sol[r_,θ_,t_,k_]:=(∑_(m=0)^k)(∑_(n=1)^k)BesselJ[m,r (Subscript[λ, m,n])]((Subscript[a, m,n])Cos[m θ]+(Subscript[b, m,n])Sin[m θ])Cos[c (Subscript[λ, m,n]) t]+(∑_(m=0)^k)(∑_(n=1)^k)BesselJ[m,r (Subscript[λ, m,n])]((Subscript[a2, m,n])Cos[m θ]+(Subscript[b2, m,n])Sin[m θ])Sin[c (Subscript[λ, m,n]) t]
solution[r_,θ_,t_]=sol[r,θ,t,5]

Plots (slow to create, but interesting):

Do[

plot1[t]=RevolutionPlot3D[solution[r,θ,t],{r,0,1},{θ,0,2Pi},PlotRange->{{-1,1},{-1,1},{-1.0,1.0}},Boxed->False,Axes->False,ImageSize->{400,Automatic}],

{t,0,5,0.125}]

Manipulate[Show[plot1[t]],{t,0,5,0.125}]