Solving the wave equation in polar coordinates

NotebookDirectory[]

Sturm-Liouville problem for θ

The SLP here is

-y''[θ]=ν y[x], y[-π]=y[π],y'[-π]=y'[π]

Here I use ν (Greek nu) for the separation of  r  and  θ  constant.  Ordinarily I would use λ, but in this case λ has already been used as a separation of time and space constant. We do know the solutions of the above SLP. However, below I show how to get the known solutions from Mathematica.

ν < 0

As usual we set ν = -(μ^2), μ>0.

soln=y[x]/.FullSimplify[DSolve[{- y''[x]+(μ^2) y[x]==0},y[x],x],μ>0][[1]]

eq1n=Collect[(soln/.{x->-π})-(soln/.{x-> π}),{C[1],C[2]}]

eq2n=Collect[(D[soln,x]/.{x->-π})-(D[soln,x]/.{x-> π}),{C[1],C[2]}]

matn={{Coefficient[eq1n,C[1]],Coefficient[eq1n,C[2]]},{Coefficient[eq2n,C[1]],Coefficient[eq2n,C[2]]}}

FullSimplify[Det[matn]]

Plot[Evaluate[FullSimplify[Det[matn]]],{μ,0,3},PlotRange->{0,5}]

Thus, no negative eigenvalues.

ν = 0

sol0=y[x]/.FullSimplify[DSolve[{- y''[x]==0},y[x],x]][[1]]

eq10=Collect[(sol0/.{x->0})-(sol0/.{x->2 π}),{C[1],C[2]}]

eq20=Collect[(D[sol0,x]/.{x->0})-(D[sol0,x]/.{x->2 π}),{C[1],C[2]}]

{{Coefficient[eq10,C[1]],Coefficient[eq10,C[2]]},{Coefficient[eq20,C[1]],Coefficient[eq20,C[2]]}}

FullSimplify[Det[{{Coefficient[eq10,C[1]],Coefficient[eq10,C[2]]},{Coefficient[eq20,C[1]],Coefficient[eq20,C[2]]}}]]

Thus, C[2]=0 and 0 is an eigenvalue with the constant function 1 as a corresponding eigenfunction.

ν > 0

As usual we set ν = (μ^2), μ>0.

solp=y[x]/.FullSimplify[DSolve[{- y''[x]-(μ^2) y[x]==0},y[x],x],μ>0][[1]]

eq1p=Collect[(solp/.{x->0})-(solp/.{x->2 π}),{C[1],C[2]}]

eq2p=Collect[(D[solp,x]/.{x->0})-(D[solp,x]/.{x->2 π}),{C[1],C[2]}]

matp={{Coefficient[eq1p,C[1]],Coefficient[eq1p,C[2]]},{Coefficient[eq2p,C[1]],Coefficient[eq2p,C[2]]}}

FullSimplify[Det[matp]]

Plot[%,{μ,0,3}]

Thus, positive eigenvalues are (m^2), m ∈{1,2,3,...}.

FullSimplify[matp/.{μ->m},m∈Integers]

Thus, both Cos[m θ]   and  Sin[m θ]  are corresponding eigenfunctions.

Sturm-Liouville problem for r

The following problem is a singular Sturm-Liouville problem

-(d/(d r))(r ((d )/(d r))y(r) )+((m^2)/r)y(r)= λ r y(r),  0<x≤1,  
   |y(r)|  bounded near r=0,    
     y(1)=0 .

The singularity comes from the fact that the coefficient (the function r) in the first term on the left hand side vanishes at r=0  and the coefficient with y(r) is not defined at 0.  This will be reflected in the solutions.  However, if we considered a drum with a hole, there would be no singularities.

Case 1: λ<0

As usual we set λ=-(μ^2), μ>0. We need Mathematica’s help in finding a solution of the differential equation -(x y'(x))'+((m^2)/x)y(x) +(μ^2)x y(x)= 0  subject to the boundary condition y(1) = 0. After we find this solution we will look for μ for which the first boundary condition y(x) bounded near x=0 is satisfied.

FullSimplify[DSolve[{-x y''[x]-y'[x]+((m^2)/x)y[x]+(μ^2) x y[x]==0,y[1]==0},y[x],x],μ>0]

Since the expression BesselY[m,-i x μ] is not defined for x=0, we conclude that none of the solutions found by Mathematica satisfy the second boundary condition. Thus, for every m=0,1,2,..., there are no negative eigenvalues.

Case 2: λ=0

Clear[y];
DSolve[{-x y''[x]-y'[x]+((m^2)/x)y[x]==0,y[1]==0},y[x],x]

It is clear that for positive m the above function is not bounded near x=0. Thus, 0 is not an eigenvalue for all positive m.

It is not clear what the above solution means for  m=0

DSolve[{-x y''[x]-y'[x]+((0^2)/x)y[x]==0,y[1]==0},y[x],x]

Since this function is not defined at x=0,  0 is not an eigenvalue when m=0.

Case 3: λ>0

As usual we set λ=(μ^2), μ>0. We first look for a fundamental set of solutions:

FullSimplify[DSolve[{-x y''[x]-y'[x]+((m^2)/x)y[x]-(μ^2) x y[x]==0},y[x],x],μ>0]

The fundamental set of solutions suggested by Mathematica is BesselJ[m,x μ]  and BesselY[m,x μ]. Let us explore these functions:

Manipulate[
Plot[
Evaluate[{BesselJ[m,x μ],BesselY[m,x μ]}/.{μ-> 1}],{x,0,20},PlotStyle->{{Thickness[0.007],Blue},{Thickness[0.007],Magenta}},PlotRange->{-4,1}
],
{{m,1},Range[0,12],ControlType->Setter}
]

Since the magenta functions are unbounded at zero they do not satisfy the first boundary condition.

Thus, the blue functions BesselJ[m, x μ] is the only solution of de which satisfies the bc at 0.  To satisfy the bs at 1 we need   BesselJ[m, 1 μ] = 0, that is

to find the positive eigenvalues we need to find μ>0  such that BesselJ[m,1 μ]=0

Manipulate[
Plot[
Evaluate[{BesselJ[m,1 μ]}],{μ,0,50},PlotStyle->{{Thickness[0.005],Blue}},Ticks->{Range[1,100,2],Automatic}
],
{{m,1},Range[0,10],ControlType->Setter}
]

These are the smallest μ-s for eigenvalues

μ/.FindRoot[BesselJ[2,μ]==0,{μ,5}]

(μ/.FindRoot[BesselJ[#,μ]==0,{μ,#+3},MaxIterations->300,AccuracyGoal->Automatic,PrecisionGoal->Infinity,WorkingPrecision->40])&/@Range[0,10]

These are good guesses for the first μ-s, but how to get the others?  By trial and error I created this table

guess=Table[Table[2.5+1.2 m+(3.1+(m/25))(n-1) ,{n,1,10}],{m,0,10}]

And this plot shows that my guesses are ok

Manipulate[
Plot[
Evaluate[{BesselJ[m,1 μ]}],{μ,0,50},PlotStyle->{{Thickness[0.005],Blue}},
Epilog->{{PointSize[0.02],Red,Point[{#,0}]&/@guess[[m+1]]}},
Ticks->{Range[1,100,2],Automatic}
],
{{m,1},Range[0,10],ControlType->Setter}
]

Now I use the guesses to get approximate solutions

evmus=Table[(μ/.FindRoot[BesselJ[k,μ]==0,{μ,#},MaxIterations->200])&/@guess[[k+1]],{k,0,10}];evmus//TableForm

Manipulate[
Plot[
Evaluate[{BesselJ[m,1 μ]}],{μ,0,50},PlotStyle->{{Thickness[0.005],Blue}},
Epilog->{{PointSize[0.02],Green,Point[{#,0}]&/@evmus[[m+1]]}},
Ticks->{Range[1,100,2],Automatic}
],
{{m,1},Range[0,10],ControlType->Setter}
]

Eigenfunctions

We now have the list of  110  μ-s in the list evmus.  The corresponding Eigenfunctions are:  
The eigenfunction corresponding to n-th  μ  corresponding to the m-th equation is

mm=3;nn=4;BesselJ[mm,evmus[[mm+1,nn]]x]

Manipulate[
Plot[
Evaluate[BesselJ[m,evmus[[m+1,n]]x]],{x,0,1},PlotStyle->{{Thickness[0.005],Blue}},PlotPoints->150,
PlotRange->{{-0.1,1+.1},{-1.1,1.1}},Ticks->{Range[0,1,.1],Automatic}
],
{{m,1},Range[0,10],ControlType->Setter},
{{n,1},Range[1,10],ControlType->Setter}
]

The eigenfunctions are mutually orthogonal:

mm=7; Chop[Table[NIntegrate[
BesselJ[mm,evmus[[mm+1,j]]x]BesselJ[mm,evmus[[mm+1,k]]x] x,{x,0,1},MaxRecursion->200,AccuracyGoal->10,PrecisionGoal->10],{j,1,10},{k,1,10}]]//MatrixForm

For later use we need the squares of the norms of the eigenfunctions:

Clear[efns]; 

efns=Chop[Table[Table[NIntegrate[
BesselJ[m,evmus[[m+1,k]]x]BesselJ[m,evmus[[m+1,k]]x] x,{x,0,1},MaxRecursion->200,AccuracyGoal->10,PrecisionGoal->10],{k,1,10}],{m,0,10}]]

Application to the unit drum

Natural modes of vibration

vp={2.3790458212203394`,-1.9562027351402624`,1.4012179843135848`}

bc=Graphics3D[{Thickness[0.006],
Line[{Cos[#],Sin[#],0}&/@Range[0,2 Pi,Pi/64]]
}, PlotRange->{{-1.1,1.1},{-1.1,1.1},{-0.07,0.07}},BoxRatios->{1,1,.4},ImageSize->300,ViewPoint->{2.3790458212203394`,-1.9562027351402624`,1.4012179843135848`}]

Manipulate[Show[ParametricPlot3D[Evaluate[
{r Cos[θ],r Sin[θ],BesselJ[m,evmus[[m+1,n]] r]*Cos[m θ]}
],{r,0,1},{θ,0,2 Pi},PlotPoints->{20,64},Mesh->False,PlotRange->{{-1.1,1.1},{-1.1,1.1},{-1.1,1.1}},Boxed->False,Axes->False,BoxRatios->{1,1,.5},ImageSize->350,ViewPoint->vp],bc],
{{m,1},Range[0,10],ControlType->Setter},{{n,1},Range[10],ControlType->Setter}]

BesselJ[0,evmus[[1,1]] 0]

mm=2;nn=1;
AbsoluteTiming[Show[ParametricPlot3D[Evaluate[N[
{r Cos[θ],r Sin[θ],BesselJ[mm,evmus[[mm+1,nn]] r]*Cos[mm θ]*Cos[evmus[[mm+1,nn]]1]}
]],{r,0,1},{θ,0,2 Pi},PlotPoints->{16,48},Mesh->False, PlotRange->{{-1.1,1.1},{-1.1,1.1},{-1,1}},Axes-> False,Boxed->False,BoxRatios->{1,1,.4},ImageSize->500,ViewPoint->vp],bc]]

5.25/0.025

Clear[ani];mm=2;nn=1;
ani=Table[Show[ParametricPlot3D[Evaluate[N[
{r Cos[θ],r Sin[θ],BesselJ[mm,evmus[[mm+1,nn]] r]*Cos[mm θ]*Cos[evmus[[mm+1,nn]]t]}
]],{r,0,1},{θ,0,2 Pi},PlotPoints->{16,48},Mesh->False, PlotRange->{{-1.1,1.1},{-1.1,1.1},{-1,1}},Axes-> False,Boxed->False,BoxRatios->{1,1,.4},ImageSize->500,ViewPoint->vp],bc],
{t,0,5.25,0.025}];

ani[[97]]

ListAnimate[ani,AnimationRate->10
]

(* Clear[mm,nn,ani];

mm=0;nn=1;

While[mm<6,
While[nn<6,
Clear[ani];
ani=Table[Show[ParametricPlot3D[Evaluate[N[
{r Cos[θ],r Sin[θ],BesselJ[mm,evmus[[mm+1,nn]] r]*Cos[mm θ]*Cos[evmus[[mm+1,nn]]t]}
]],{r,0,1},{θ,0,2 Pi},PlotPoints->{16,48},Mesh->False, PlotRange->{{-1.1,1.1},{-1.1,1.1},{-1,1}},Axes-> False,Boxed->False,BoxRatios->{1,1,.4},PlotLabel->TableForm[{"t=",NumberForm[t,{3,3}]},TableDirections->Row,TableSpacing->0.2],ImageSize->500,ViewPoint->vp],bc],
{t,0,5.25,0.025}];Export[StringJoin["DrumNMV",ToString[mm],ToString[nn],"s1.gif"],ani[[1]],"GIF",{"ImageSize"->500}];
Export[StringJoin["DrumNMV",ToString[mm],ToString[nn],"ani.gif"],ani,"GIF","AnimationRepetitions"->0,"ImageSize"->500,"DisplayDurations"->0.15] ;
nn++];nn=1;mm++] *)

(* 

SetDirectory["C:\\Dropbox\\Work\\myweb\\Courses\\Math_pages\\Math_430"]; 

Export["DrumNMV21s1.gif",ani[[1]],"GIF",{"ImageSize"->800}];

Export["DrumNMV21ani.gif",ani,"GIF","AnimationRepetitions"->0,"ImageSize"->800,"DisplayDurations"->0.1] 

*)

Vibrations initiated by initial velocity

Clear[gr,gt];gr[r_]=2(1-r)((E^(-100(((r-.5))^2))));gt[θ_]=Sin[(E^(-50 (((θ-π))^2)))*(π/2)]+.0000000000000000000000000001

gr[1]

Show[ParametricPlot3D[Evaluate[ 
{r Cos[θ],r Sin[θ],gr[r]gt[θ]}
 ],{r,0,1},{θ,0, 2Pi},PlotPoints->{10,28},Mesh->False, PlotRange->{{-1.1,1.1},{-1.1,1.1},{-1,1}},BoxRatios->{1,1,1},ImageSize->350,ViewPoint->Dynamic[vp]],bc]

efns[[1]]

Clear[coefr1];

coefr1=Table[(1/(efns[[1,n]]))NIntegrate[gr[r] BesselJ[0,evmus[[1,n]] r]r,{r,0,1},Method->{Automatic,"SymbolicProcessing"->0}],{n,1,10}]

Clear[coefr];

coefr=Table[Table[(1/(efns[[m+1,n]]))NIntegrate[gr[r] BesselJ[m,evmus[[m+1,n]] r]r,{r,0,1},Method->{Automatic,"SymbolicProcessing"->0}],{n,1,10}],{m,0,10}]

Clear[coefct];

coefct=Table[(1/Pi)NIntegrate[gt[t] Cos[m t],{t,0,2 Pi},Method->{Automatic,"SymbolicProcessing"->0}],{m,1,10}]

coef1t=(1/(2Pi))NIntegrate[gt[t]  1,{t,0,2 Pi},Method->{Automatic,"SymbolicProcessing"->0}]

Clear[coefst];

coefst=Chop[Table[(1/Pi)NIntegrate[gt[t] Sin[m t],{t,0,2 Pi},MaxRecursion->200,
PrecisionGoal->16,WorkingPrecision->10,AccuracyGoal->12],{m,1,10}]]

Clear[approx];
approx[r_,θ_]=coef1t((∑_(n=1)^10)coefr[[1,n]] BesselJ[0,evmus[[1,n]] r])+(∑_(m=1)^10)coefct[[m]]Cos[m θ]((∑_(n=1)^10)coefr[[m+1,n]] BesselJ[m,evmus[[m+1,n]] r]);

approx[.3,.4]

This is the approximation of the initial velocity

Show[ParametricPlot3D[Evaluate[N[
{r Cos[θ],r Sin[θ],approx[r,θ]}
]],{r,0,1},{θ,0,2 Pi},PlotPoints->{10,32},Mesh->False, PlotRange->{{-1.1,1.1},{-1.1,1.1},{-1,1}},BoxRatios->{1,1,1},ImageSize->350,ViewPoint->Dynamic[vp]],bc]

And this is the difference between the given velocity and our approximation

ParametricPlot3D[Evaluate[N[
{r Cos[θ],r Sin[θ],approx[r,θ]-gr[r]gt[θ]}
]],{r,0,1},{θ,0,2 Pi},PlotPoints->{10,32},Mesh->False, PlotRange->{{-1.1,1.1},{-1.1,1.1},{-.1,.1}},BoxRatios->{1,1,1},ImageSize->350,ViewPoint->vp]

Next we define our solution:

Clear[uu];
uu[r_,θ_,t_]=N[coef1t((∑_(n=1)^10)((coefr[[1,n]])/(evmus[[1,n]])) BesselJ[0,evmus[[1,n]] r]Sin[evmus[[1,n]]t])+(∑_(m=1)^10)coefct[[m]]Cos[m θ]((∑_(n=1)^10)((coefr[[m+1,n]])/(evmus[[m+1,n]])) BesselJ[m,evmus[[m+1,n]] r]Sin[evmus[[m+1,n]] t])];

uu[.3,2.3,2]

AbsoluteTiming[Show[ParametricPlot3D[Evaluate[N[
{r Cos[θ],r Sin[θ],uu[r,θ,.41]}
]],{r,0,1},{θ,0,2 Pi},PlotPoints->{16,48},Mesh->False, PlotRange->{{-1.1,1.1},{-1.1,1.1},{-0.07,0.07}},BoxRatios->{1,1,.4},Boxed->False,Axes->False,PlotLabel->TableForm[{"t=",NumberForm[t,{3,2}]},TableDirections->Row,TableSpacing->0.2],ImageSize->400,ViewPoint->{2.3790458212203394`,-1.9562027351402624`,1.4012179843135848`}],bc]]

(* 

This takes long to evaluate. 

Clear[aniuu];
aniuu=Table[Show[ParametricPlot3D[Evaluate[N[
{r Cos[θ],r Sin[θ],uu[r,θ,t]}
]],{r,0,1},{θ,0,2 Pi},PlotPoints->{20,64},Mesh->False, PlotRange->{{-1.1,1.1},{-1.1,1.1},{-0.07,0.07}},BoxRatios->{1,1,.4},Boxed->False,Axes->False,PlotLabel->TableForm[{"t=",NumberForm[t,{3,2}]},TableDirections->Row,TableSpacing->0.2],ImageSize->800,ViewPoint->{2.3790458212203394`,-1.9562027351402624`,1.4012179843135848`}],bc],
{t,0.,7.,0.05}]; *)

(* Show[aniuu[[5]],ImageSize->300] *)

(* ListAnimate[aniuu,AnimationRate->8
]  *)

(*  SetDirectory["C:\\Dropbox\\Work\\myweb\\Courses\\Math_pages\\Math_430"]; 

Export["DrumVs1.gif",aniuu[[1]],"GIF",{"ImageSize"->800}];

Export["DrumVani.gif",aniuu,"GIF","AnimationRepetitions"->0,"ImageSize"->800,"DisplayDurations"->0.3] 

*)