Partial Differential EquationsHigher-dimensional PDE: Vibrating circular membranes and Bessel functions.Anton DzhamayDepartment of MathematicsThe University of MichiganAnn Arbor, MI 48109wPage: http://www.math.lsa.umich.edu/~adzhamemail: adzham@umich.eduCopyright \251 2004 by Anton DzhamayAll rights reserved

In this worksheet we consider some examples of vibrating circular membranes. Such membranes are described by the two-dimensional wave equation. Circular geometry requires the use of polar coordinates, which in turn leads to the Bessel ODE, and so the basic solutions obtained by the method of separations of variables (product solutions or standing waves) are described with the help of Bessel functions.

Some packages that we use in this worksheet:restart: with(plottools): with(plots):

Graphs of a few Bessel functions of the first and second kind:display([seq( display(plot([BesselJ(m,x),BesselY(m,x)],x=0..20,thickness=3, color=[red,blue],view=-1..1, title=sprintf("Bessel functions J_%d and Y_%d",m,m))),m=0..5)], insequence=true);For large values of NiMlInpH Bessel functions can be approximated by the usual trigonometric functions with decreasing amplitude:BesselJa:=(m,z)->sqrt(2/(Pi*z))*cos(z-Pi/4-m*Pi/2): BesselYa:=(m,z)->sqrt(2/(Pi*z))*sin(z-Pi/4-m*Pi/2):display([seq( display(plot([BesselJ(m,x),BesselJa(m,x)],x=0..20,thickness=3, color=[red,blue],view=-1..1, title=sprintf("Bessel function J_%d and its approximation",m))),m=0..5)], insequence=true);display([seq( display(plot([BesselY(m,x),BesselYa(m,x)],x=0..20,thickness=3, color=[red,blue],view=-1..1, title=sprintf("Bessel function Y_%d and its approximation",m))),m=0..5)], insequence=true);

Let us introduce a short-hand notation for the NiMlIm5H-th zero of the Bessel function NiMmJSJKRzYjJSJtRw==:z:=BesselJZeros:z(m,n);First we define the spatial eigenfunction NiMtJiUkUGhpRzYjNiQlIm1HJSJuRzYkJSJyRyUmdGhldGFH of the Dirichlet boundary problem:PhiC:=unapply(BesselJ(m,z(m,n)*r/a)*cos(m*theta),m,n):'Phi[c][m,n](x,y)'=PhiC(m,n);PhiS:=unapply(BesselJ(m,z(m,n)*r/a)*sin(m*theta),m,n):'Phi[s][m,n](r,theta)'=PhiS(m,n);The corresponding eigenvalue NiMmJSdsYW1iZGFHNiM2JCUibUclIm5H is lambda:=unapply((z(m,n)/a)^2,m,n):'lambda[m,n]'=lambda(m,n);In this worksheet we consider only initial displacements, and thereforeT:=unapply(cos(c*sqrt(lambda(m,n))*t),m,n):T(m,n);and the product solutions areuC:=unapply(PhiC(m,n)*(T(m,n)),m,n):'u[c][m,n](r,theta,t)'=uC(m,n); uS:=unapply(PhiS(m,n)*(T(m,n)),m,n):'u[s][m,n](r,theta,t)'=uS(m,n);The period of (NiQlIm1HJSJuRw==)-oscillation isP:=unapply( (2*Pi)/(c*sqrt(lambda(m,n))),m,n):P(m,n);We also introduce the following procedure that will allow us to plot solutions together with nodal curves:addcoords(z_cylindrical,[z,r,theta],[r*cos(theta),r*sin(theta),z]); wave:=proc(sol) display([ contourplot3d(sol,r=0..1.01*a,theta=0..2*Pi,contours=[0],color=red, coords=z_cylindrical,thickness=3,numpoints=600), plot3d(sol,r=0..1.01*a,theta=0..2*Pi,coords=z_cylindrical,shading=zhue)],scaling=constrained,axes=boxed); end:

c:=1:a:=2:The "lowest energy" mode has the eigenvalue'lambda(0,1)'=lambda(0,1);and the corresponding eigenfunction is obtained by mapping the first zero of NiMmJSJKRzYjIiIh to the radius of the circle: 'PhiC[0,1]'=PhiC(0,1);animate(wave,[uC(0,1)],t=0..P(0,1));Here we map the second zero NiMmJSJ6RzYkIiIhIiIj of NiMmJSJKRzYjIiIh to the radius NiMlImFH of the circle:animate(wave,[uC(0,2)],t=0..P(0,2));And the third:animate(wave,[uC(0,3)],t=0..P(0,3));Taking NiMvJSJtRyIiIg== introduces the cosine and sine terms:'PhiC[1,1]'=PhiC(1,1);'PhiS[1,1]'=PhiS(1,1);animate(wave,[uC(1,1)],t=0..P(1,1));animate(wave,[uS(1,1)],t=0..P(1,1));animate(wave,[uC(1,1)+uS(1,1)],t=0..P(1,1));One can easily identify the type of the product solution from the picture of the nodal curve: three radial lines below tell us that NiMvJSJtRyIiJA==, two circular lines tell us that NiMvJSJuRyIiIw==, and since one of the nodal lines coinsides with the NiMlInhH-axis (given by NiMvJSZ0aGV0YUciIiE=), we know that the NiMlJnRoZXRhRw==-term is a sine function:display(animate(wave,[uS(3,2)],t=0..P(3,2)),insequence=true,orientation=[-90,0]);Choosing eigenfunctions that have different eigenvalues results in "moving waves":animate(wave,[uS(1,1)+0.2*uC(0,3)],t=0..2*P(1,1)*P(0,3));animate(wave,[uS(2,2)+0.3*uC(1,4)],t=0..2*P(2,2)*P(1,4));animate(wave,[uC(0,1)+0.8*uS(1,2)+0.3*uC(1,1)],t=0..P(0,1)*P(1,1)*P(1,2));Clean-up:a:='a':c:='c':

Walter A. Strauss, Partial Differential Equations: An Introduction, Wiley, 1992Richard Haberman, Elementary Applied Partial Differential Equations, 3rd edition, Prentice HallStanley J. Farlow, Partial Differential Equations for Scientists and Engineers, Dover 1982

"While every effort has been made to validate the solutions in this worksheet, Waterloo Maple Inc. and the contributors are not responsible for any errors contained and are not liable for any damages resulting from the use of this material."