A Vibrating Circular Membrane

Prepared by Dennis W. Brewer

Introduction

The wave equation in polar coordinates is (Subscript[u, tt])=(a^2)((Subscript[u, rr])+(1/r)(Subscript[u, r])+(1/(r^2))(Subscript[u, θθ])).  If we assume radial symmetry of the initial conditions then u  is independent of θ  and the problem reduces to (Subscript[u, tt])=(a^2)((Subscript[u, rr])+(1/r)(Subscript[u, r])) with boundary condition u(c,t)=0,  t>0 and initial conditions u(r,0)=f(r), (Subscript[u, t])(r,0)=g(r),  0≤r<c.  
Proceeding by the product method, we seek a solution of the form u(r,t)=R(r)T(t).  Substituting into the PDE and separating the variables yields ((T'')/((a^2)T))=((R'')/R)+((R')/rR)=-λ which gives to ordinary differential equations T''+(a^2)λT=0 and rR''+R'+λrR=0.  Since we expect the solutions to oscillate in time, we take λ=(α^2)>0 so that T(t)=(Subscript[c, 1])cos(aαt)+(Subscript[c, 2])sin(aαt).  
Now the equation rR''+R'+(α^2)rR=0 is related to Bessel's equation of order zero.  Here is the general solution using Mathematica.

DSolve[r*R''[r]+R'[r]+α^2*r*R[r]==0,R[r],r]

The functions (Subscript[J, 0]) and (Subscript[Y, 0]) are the Bessel functions of order zero of the first and second kind, respectively.  Here are plots of these functions.

Plot[{BesselJ[0,r],BesselY[0,r]},{r,0,10}]

The function (Subscript[Y, 0])(r) approaches negative infinity as r→(SuperPlus[0]) so we will take the solution R(r)=(Subscript[c, 3])(Subscript[J, 0])(αr).   Now from the boundary condition on the circumference of the circle we must have R(c)=0 so that (Subscript[J, 0])(αc)=0.  From the plot above, it appears that (Subscript[J, 0]) has infinitely many positive zeros.  We will denote these zeros by (Subscript[j, 1]),(Subscript[j, 2]),(Subscript[j, 3]), ....  In Mathematica the nth zero of (Subscript[J, 0])is denoted by BesselJZero[0,n].    Here is a list of the first ten zeros.

Table[N[BesselJZero[0,n]],{n,1,10}]

Returning to our solution, the eigenvalues are λ=((((Subscript[α, n])))^2) where (Subscript[α, n])=((Subscript[j, n])/c), n=1,2,3,....  Applying superposition to the product solutions yields a final solution of the form

({{({{u(r,t)=(∑_(n=1)^∞)((Subscript[A, n])cos (((Subscript[aj, n])t)/c)+(Subscript[B, n])sin (((Subscript[aj, n])t)/c))(Subscript[J, 0])((((Subscript[j, n])r)/c))}})}})

Note that the fundamental period is ((2πc)/(Subscript[aj, 1])) and the fundamental pitch is (Subscript[f, 0])=((Subscript[aj, 1])/(2πc)).  For example, if the radius of the drum, c, is doubled, then the fundamental pitch is half as much, or one octave lower than before.  The first overtone is at a frequency of (Subscript[f, 1])=((Subscript[aj, 2])/(2πc)) so that the ratio of the two pitches is ((Subscript[j, 2])/(Subscript[j, 1])).  Using Mathematica, the ratios of the overtones to the fundamental pitch are

Table[N[BesselJZero[0,n]/BesselJZero[0,1]],{n,1,10}]

Note that unlike the vibrating string, these are not integers!  Thus vibrating circular membranes sound different from vibrating strings!  Suppose the fundamental pitch is tuned to a "C" at 512 Hertz.  Then the frequencies of the overtones are

Table [512*N[BesselJZero[0,k]/BesselJZero[0,1]],{k,2,8}]

The first overtone lies in the crack between a D and E flat on the piano.  The second overtone is close to a B flat.  If you have a speaker, you can play the fundamental pitch and the first and second overtones together.

Play[Sin[512*2*Pi*t]+Sin[1175.25*2*Pi*t]+Sin[1842.42*2*Pi*t],{t,0,5}]

Compare this chord with the sound of a vibrating string on the same fundamental pitch.

Play[Sin[512*2*Pi*t]+Sin[1024*2*Pi*t]+Sin[1536*2*Pi*t],{t,0,5}]

Next we consider the effect of the initial conditions which which the membrane is set in motion.  From the initial conditions we obtain the formulas

({{({{u(r,0)=(∑_(n=1)^∞)(Subscript[A, n])(Subscript[J, 0])((((Subscript[j, n])r)/c))}})=f(r), 0≤r<c
({{(Subscript[u, t])(r,0)=(∑_(n=1)^∞)((Subscript[ax, n])/c)(Subscript[B, n])(Subscript[J, 0])((((Subscript[j, n])r)/c))}})=g(r), 0≤r<c}})

These are called Fourier-Bessel Series.   We can compute the cofficients of a Fourier-Bessel series by taking advantage of the orthogonality of the eigenfunctions (Subscript[J, 0])((((Subscript[j, n])r)/c)) on the interval 0<r<c with respect to the weight function p(r)=r.    The following table shows the orthogonality with c=1.

Table[Integrate[r*BesselJ[0,BesselJZero[0,m]*r]*BesselJ[0,BesselJZero[0,n]*r],{r,0,1}],{m,1,3},{n,1,3}]

Mathematica computes a general formula for the inner product of an eigenfunction with itself, namely,

Integrate[r*BesselJ[0,BesselJZero[0,n]*r/radius]^2,{r,0,radius}]

Since (Subscript[J, 0])((Subscript[j, n]))=0, this simplifies to (∫_0^c)((r((Subscript[J, 0])((((Subscript[j, n])r)/c))))^2)dr=(1/2)(((c^2)((Subscript[J, 1])((Subscript[j, n]))))^2)where (Subscript[J, 1]) denotes the Bessel function of order one, denoted by BesselJ[1,x] in Mathematica.    Therefore, using the orthogonality of the eigenfunctions we arrive at the following formulas for the coefficients of the Fourier-Bessel series.

({{({{(Subscript[A, n])=(2/((c^2)(Subscript[J, 1]^2)((Subscript[j, n]))))(∫_0^c)(Subscript[rJ, 0])((((Subscript[j, n])r)/c))f(r)dr
(Subscript[B, n])=(2/((Subscript[ax, n])(Subscript[cJ, 1]^2)((Subscript[j, n]))))(∫_0^c)(Subscript[rJ, 0])((((Subscript[j, n])r)/c))g(r)dr}})}})

Example  (Subscript[u, tt])=(a^2)((Subscript[u, rr])+(1/r)(Subscript[u, r])), 0<r<c,t>0, u(c,t)=0 , t>0, u(r,0)=0, (Subscript[u, t])(r,0)=g(r), 0≤r<1, where g(r)=-(Subscript[v, 0]) for 0≤r<b and g(r)=0 for b≤r<1.

The motivation behind these initial conditions is that the membrane is set in motion from its equiliibrium position by striking it in the center.
We will use the following parameter values and set the number of terms in the Fourier-Bessel series.

a=1;c=1;v0=1;b=.1;ncoeffs=5;

Here are the Fourier-Bessel coefficients corresponding to the initial velocity function in the example.

Bcoeffs=Table[(2/(a*BesselJZero[1,n]*c*BesselJ[1,BesselJZero[0,n]]^2))*NIntegrate[r*(-v0)*BesselJ[0,BesselJZero[0,n]*r/c],{r,0,b}],{n,1,ncoeffs}]

Let's make a numerical table of the zeros of the Bessel function to speed up the computations later on.

j=Table[N[BesselJZero[0,n]],{n,1,ncoeffs}]

Here is the resulting approximation of the initial velocity.

Plot[Bcoeffs.Table[(a*j[[n]]/c)*BesselJ[0,j[[n]]*r/c],{n,1,ncoeffs}],{r,0,c}, PlotRange->{-1.2,.2}]

Here is a plot of the solution along a fixed line from the center of the circle as a function of time.

Manipulate[Plot[Bcoeffs.Table[Sin[a*j[[n]]*t/c]*BesselJ[0,j[[n]]*r/c],{n,1,ncoeffs}],{r,0,c},PlotRange->.07],{t,0,2*Pi*c/(a*j[[1]])}]

Here is a three-dimensional plot of the vibrating circular membrane.  Step through time slowly to give the computer time to catch up!

Manipulate[RevolutionPlot3D[Bcoeffs.Table[Sin[a*j[[n]]*t/c]*BesselJ[0,j[[n]]*r/c],{n,1,ncoeffs}],{r,0,c},{phi,0,2*Pi}, PlotRange->{All,All,{-.05,.05}}],{t,0,2*Pi*c/(a*j[[1]])}]

Here are plots of several standing waves that can make up the full vibration of the circular membrane.

Manipulate[RevolutionPlot3D[Bcoeffs[[1]]
*Sin[a*j[[1]]*t/c]*BesselJ[0,j[[1]]*r/c],{r,0,c},{phi,0,2*Pi}, PlotRange->{All,All,{-.01,.01}}],{t,0,2*2*Pi*c/(a*j[[1]])}]

Manipulate[RevolutionPlot3D[Bcoeffs[[2]]
*Sin[a*j[[2]]*t/c]*BesselJ[0,j[[2]]*r/c],{r,0,c},{phi,0,2*Pi}, PlotRange->{All,All,{-.015,.015}}],{t,0,2*2*Pi*c/(a*j[[1]])}]

Manipulate[RevolutionPlot3D[Bcoeffs[[3]]
*Sin[a*j[[3]]*t/c]*BesselJ[0,j[[3]]*r/c],{r,0,c},{phi,0,2*Pi}, PlotRange->{All,All,{-.015,.015}}],{t,0,2*2*Pi*c/(a*j[[1]])}]

Here is the sound of the vibrating circular membrane with the frequencies scaled up to the normal range of human hearing.

Play[Bcoeffs.Table[Sin[a*j[[n]]*500*t/c]*BesselJ[0,j[[n]]*0/c],{n,1,ncoeffs}],{t,0,3}]