Circular Membrane with 
Fixed Boundary (Simplified)

Statement of the problem

In this file I will consider the following wave equation 
  
     (((∂^2)u)/(∂(t^2)))=(c^2) ((((∂^2)u)/(∂(x^2)))+(((∂^2)u)/(∂(y^2)))),  (x,y) ∈ D,   0 ≤ t < ∞.  
  
Here D is the unit cirlce in x y-plane (or, rather, r θ-plane): 

  D = {(x,y) ∈ (ℝ^2) |  (x^2)+(y^2)  ≤ 1}={(r,θ) ∈ (ℝ^2) |  r  ≤ 1,  0⩽θ⩽ 2π} 
 
and  u(x,y,t)  is the function describing the shape of the membrane at time  t .  The equation is subject to the boundary condition  u(x,y,0)= 0   whenever  (x^2)+(y^(2 ))= 1.

Clearly it is easier to deal with this equation in polar coordinates.

Let 
 v(r,θ,t) = u(r Cos[θ],r Sin[θ],t).   
Then 
((∂v)/(∂r))=((∂u)/(∂x))Cos[θ]+((∂u)/(∂y))Sin[θ]    and  ((∂v)/(∂θ))=-((∂u)/(∂x))r Sin[θ]+((∂u)/(∂y))r Cos[θ] .  
Solving for  ((∂u)/(∂x)) and  ((∂u)/(∂y))  we get 
((∂u)/(∂x))=((∂v)/(∂r))  Cos[θ]-(1/r)((∂v)/(∂θ))Sin[θ]  and  ((∂u)/(∂y))=((∂v)/(∂r)) Sin[θ]+(1/r)((∂v)/(∂θ)) Cos[θ] . 
Taking the second partial derivatives of   v(r,θ,t) we get  
(((∂^2)v)/(∂(r^2)))=(((∂^2)u)/(∂(x^2)))(Cos^2)[θ]+(((∂^2)u)/(∂y∂x))Sin[2θ]+(((∂^2)u)/(∂(y^2)))(Sin^2)[θ]  
and  
(((∂^2)v)/(∂(θ^2)))=(((∂^2)u)/(∂(x^2)))(r^2) (Sin^2)[θ]-(((∂^2)u)/(∂y∂x))(r^2) Sin[2θ] + (((∂^2)u)/(∂(y^2)))(r^2)(Cos^2)[θ]-((∂u)/(∂x))r Cos[θ] -((∂u)/(∂y))r Sin[θ] .  
The last equation simplifies to 

 (1/(r^2))((((∂^2)v)/(∂(θ^2)))+r ((∂v)/(∂r)))=(((∂^2)u)/(∂(x^2)))(Sin^2)[θ]-(((∂^2)u)/(∂y∂x))Sin[2θ] + (((∂^2)u)/(∂(y^2)))(Cos^2)[θ].  

Adding the last equation and the equation 
 (((∂^2)v)/(∂(r^2)))=(((∂^2)u)/(∂(x^2)))(Cos^2)[θ]+(((∂^2)u)/(∂y∂x))Sin[2θ]+(((∂^2)u)/(∂(y^2)))(Sin^2)[θ] ,
we get 

(((∂^2)u)/(∂(x^2)))+(((∂^2)u)/(∂(y^2)))= (1/(r^2))(((∂^2)v)/(∂(θ^2)))+(((∂^2)v)/(∂(r^2))) +(1/r)((∂v)/(∂r)) .

Thus, the two dimensional wave equation in polar coordinates becomes 

    (((∂^2)v)/(∂(t^2)))=(c^2) ((1/(r^2))(((∂^2)v)/(∂(θ^2)))+(((∂^2)v)/(∂(r^2))) +(1/r)((∂v)/(∂r))),  
     0 ≤ r ≤ 1, 0 ≤ θ <2π,   0 ≤ t < ∞,
	v(1,θ,t) = 0,  0 ≤ θ <2π,   0 ≤ t < ∞ , 
	v(r,θ,0) = f(r,θ), ((∂v)/(∂t))(r,θ,0)=g(r,θ).

The separation of variables

As we discussed in class to solve this problem by separation of variables  we look for the solution of the form    v(r,θ,t) = A(r,θ) B(t).  The partial differential equation becomes: 

A(r,θ) B''(t) = (c^2) ((1/(r^2))(Subscript[A, θθ])(r,θ) +(Subscript[A, r r])(r,θ) + (1/r)(Subscript[A, r])(r,θ))B(t),
   0 ≤ r ≤ 1, 0 ≤ θ <2π,   0 ≤ t < ∞,

or equivalently,

((B''(t))/((c^2)B(t)))=(((1/(r^2))(Subscript[A, θθ])(r,θ) +(Subscript[A, r r])(r,θ) + (1/r)(Subscript[A, r])(r,θ))/(A(r,θ)))=-λ, 
       0 < r ≤ 1, 0 ≤ θ <2π,   0 ≤ t < ∞.                  (1)

The first set of boundary conditions leads to 

A(1,θ) = 0    and a natural condition that 

A(r,0) = A(r,2π)   and  (Subscript[A, θ])(r,0) = (Subscript[A, θ])(r,2π). 

Since the independent variables in (1) are separated we have the following equations:  For the unknown function A :

(1/(r^2))(Subscript[A, θθ])(r,θ) +(Subscript[A, r r])(r,θ) + (1/r)(Subscript[A, r])(r,θ)=-λ A(r,θ),  
               0 < r ≤ 1, 0 ≤ θ ≤ 2π,   
             A(1,θ) = 0 
            A(r,0) = A(r,2π),  
            (Subscript[A, θ])(r,0) = (Subscript[A, θ])(r,2π).       
                                                                                                                        (2)

The above boxed problem is called a boundary eigenvalue problem.  
For the unknown function B we get:

-B''(t) = λ (c^2) B(t),   0 ≤ t < ∞ .

Note that the unknown quantity   λ  in the equation for  B  is green since we have to solve the boundary eigenvalue problem (2) for both A  and   λ .

A simplification

In this notebook we will assume that the drum is radialy symmetric. That is the function A(r,θ) does not depend on θ. Thus, instead of A(r,θ), we can simply write A(r). In this case the problem (2) simplifies to:

(Subscript[A, r r])(r) + (1/r)(Subscript[A, r])(r)=-λ A(r),    0 < r ≤ 1,    
            A(1) = 0     and   
          |A(r)|  < +∞ ,    0 < r ≤ 1 .
                                                                                                                                (2s)

This is a singular Sturm-Liouville eigenvalue problem, with singularity at 0.

Note that this simplification implies that the initial shape of the drum and the initial velocity given to the drum can not depend on θ. Thus the given functions f and g  depend only on r.

Solving the boundary eigenvalue problem for the function A and the eigenvalue parameter  λ :

A'(r) + (1/r)A(r)=-λ A(r),  0 < r ≤ 1,    
            A(1) = 0   and   
         |A(r)|   < +∞,    0 < r ≤ 1 .
                                                                                                                                (2s)

Lets write this problem in the standard Sturm-Liouville form:

-(r A'(r))'=λ r A(r),  0 < r ≤ 1,    
            A(1) = 0   and   
          |A(r)|< +∞ ,    0 < r ≤ 1.
                                                                                                                                (2sl)

Now we apply the standard procedure to find eigenvalues and eigenfunctions.  In fact we solved this problem in our notebook StuLio4a.nb.  Lets copy the solution from that notebook as our conclusion.

Expand[DSolve[r A''[r]+A'[r]+(μ^2)r A[r]==0,A[r],r][[1]][[1]][[2]]]

Read StuLio4a.nb to understand the conclusion below.

Conclusion :

We have found the squareroots of the eigenvalues of the problem (2s) by solving the following equation:  

BesselI[0,I μ] = 0

with μ > 0.  Looking at the graph:

Plot[BesselI[0,I μ],{μ,0,35},PlotRange->{-1.1,1.1}]

Lets find the approximations:

guess={2,6,9,12,15,18,21,24,27,31};
tevs=N[Chop[Table[FindRoot[BesselI[0,I μ]==0,{μ,guess[[n]]},MaxIterations->30,WorkingPrecision->20][[1]][[2]],{n,1,10}]],12]

Plot[{BesselI[0,I μ]},{μ,0,35},PlotRange->{-1,1},PlotStyle->{{RGBColor[1,0,0]}},
Ticks->{
			Table[{N[tevs[[n]],3],N[tevs[[n]],3],.3},{n,1,10}],Automatic}]

The eigenfunctions are

A[r_,n_]:=BesselI[0,I r μ]/.μ->tevs[[n]]

Table[A[r,n],{n,1,10}]

Manipulate[Plot[A[r,n],{r,0,1},PlotRange->{-1,1}],{n,1,10,1}]

Check orthogonality:

Off[NIntegrate::precw];
Chop[Table[NIntegrate[Expand[A[r,n]* A[r,m]*r],{r,0,1},MaxRecursion->20,
PrecisionGoal->12,WorkingPrecision->20,AccuracyGoal->12],{n,1,2},{m,1,4}]]//TableForm

Off[NIntegrate::precw];
ns=Chop[Table[NIntegrate[Expand[A[r,n]* A[r,n]*r],{r,0,1},MaxRecursion->20,
PrecisionGoal->12,WorkingPrecision->20,AccuracyGoal->12],{n,1,10}],(10^(-9))]

Solving the equation for the function B:

In the equation for  B  the quantity λ=(μ^2), and we found the first ten μ in the Concusion section.  They are in the table tevs. To find  B  we have to solve:

-B''(t) = (μ^2) (c^2) B(t),   0 ≤ t < ∞ .

Thus, the solutions that we find for  B   that correspond to the eigenvalue are
        Cos[c μ t]  and   Sin[c μ t].

Define:

c=1

(Subscript[B, p])[t_,n_]:=Cos[c μ t]/.μ->tevs[[n]];
(Subscript[B, v])[t_,n_]:=Sin[c μ t]/.μ->tevs[[n]];

Few solutions of the problem:

Finally, we summerize our investigations by concluding that we have found "few" solutions:  for each μ in the table (which can be as large as we want), we found one solutions for A(r) and two solutions for B(t): 

A[r,n]*(Subscript[B, p])[t,n]     and     A[r,n]*(Subscript[B, v])[t,n] ,    n = 1, 2,... .

It is interesting to write the exact formulas for the "few" solutions: 

BesselI[0,I r μ] *Cos[c μ t]   and  BesselI[0,I r μ] *Sin[c μ t]   

where   μ  is the solution of the equation  BesselI[0,I μ] = 0.

Here n  relates to the eigenvalue in the table of the squareroots of eigenvalues.  To plot these functions we have to use parametric plots of surfaces. 
Lets explore these functions.

n=3;c=1;
ParametricPlot3D[{r Cos[θ],r Sin[θ],A[r,n]*(Subscript[B, p])[0,n]},{r,0,1},{θ,0,2π},
		PlotRange->{{-1,1},{-1,1},{-1.2,1.2}},PlotPoints->{30,30},BoxRatios->{1,1,.6}]

n=3;c=1;
Manipulate[
	ParametricPlot3D[{r Cos[θ],r Sin[θ],A[r,n]*(Subscript[B, p])[t,n]},{r,0,1},{θ,0,2π},
			PlotRange->{{-1,1},{-1,1},{-1.2,1.2}},PlotPoints->{30,30},BoxRatios->{1,1,.6}],{t,0,((2π)/(c tevs[[n]])),((2π*(1/25))/(c tevs[[n]]))}]

From Few Many (FFM) and the solution of the problem:

From these "few" solutions the superposition principle allows us to create "many" solutions:  

v(r,θ,t)  = (∑_(n=1)^∞)(Subscript[a, n ])A[r,n]* (Subscript[B, p])[t,n ] + (∑_(n=1)^∞)(Subscript[b, n ]) A[r,n]*(Subscript[B, v])[t,n] 

where (Subscript[a, n]) and (Subscript[b, n ]) are arbitrary constants.

We can determine  (Subscript[a, n]) and (Subscript[b, n ]) such that 

f(r) = v(r,θ,0) =(∑_(n=1)^∞)(Subscript[a, n  ]) A[r,n] 

and 

g(r) = (∂/(∂t))v(r,θ,0) =(∑_(n=1)^∞)(Subscript[b, n ])c (Subscript[μ, n])A[r,n] .

Using the orthogonality of   A[r,n] , n = 1,2, ... , on  0 ⩽  r ⩽  1, we calculate, for n=1,2, ... ,  : 

(Subscript[a, n ])=      (((∫_0^1)f(r) A[r,n]  rdr)/((∫_0^1) A[r,n]  A[r,n] r dr))

and   

(Subscript[b, n]) =   (1/(c (Subscript[μ, n])))(∫_0^1)((g(r) A[r,n]  rdr)/((∫_0^1) A[r,n]  A[r,n] r dr)).

Mathematica implementation of the solution:

g[r_]=((E^(-100(r^2)))-(E^(-100)))

Plot[g[r],{r,-1,1},PlotRange->{0,1}]

(Subscript[a, n_])[f_]:=(1/(ns[[n]]))(∫_0^1)N[f[r] A[r,n]] r \[DifferentialD]r;

(Subscript[b, n_])[g_]:=(1/(ns[[n]])) *(1/(c tevs[[n]]))(∫_0^1)g[r] A[r,n]r  \[DifferentialD]r;

Off[NIntegrate::precw];
Tb=Chop[Table[(1/(ns[[n]])) *(1/(c tevs[[n]]))NIntegrate[Expand[g[r]* A[r,n]*r],{r,0,1},MaxRecursion->20,
PrecisionGoal->8,WorkingPrecision->20,AccuracyGoal->8],{n,1,10}],(10^(-9))];

M=10;c=1;
v[r_,t_]:=Chop[(∑_(n=1)^M)Tb[[n]]*A[r,n]*(Subscript[B, v])[t,n]]

Plot3D[v[r,t],{r,0,1},{t,0,3}]

Plot[{g[r],Evaluate[(D[v[r,t],t]/.t->0)]},{r,0,1},PlotRange->{0,1}]

Manipulate[
	ParametricPlot3D[{r Cos[θ],r Sin[θ],v[r,t]},{r,0,1},{θ,0,2π},
			PlotRange->{{-1,1},{-1,1},{-.05,.05}},PlotPoints->{25,25},BoxRatios->{1,1,.6}],{t,0,((2π)/2.40),((2π)/2.40)*(1/20)}]