Solution of wave equation on a disk:

We consider the general wave equation on a disk of radius R
	(Subscript[u, tt]) = (c^2) (((Subscript[(r (Subscript[u, r])), r])/r) +((Subscript[u, θθ])/(r^2)))
subject to the boundary condition
	u(R,θ,t)=0
and initial conditions
	u(r,θ,0)=f(r,θ),
	(Subscript[u, t])(r,θ,0)=g(r,θ).
If the problem is circularly symmetric then the PDE simplifies to
	(Subscript[u, tt]) = (c^2) ((Subscript[(r (Subscript[u, r])), r])/r).
If, in addition, we assume that the initial velocity is zero, then the boundary and initial conditions become
	u(R,t)=0,
	u(r,0)=f(r),
	(Subscript[u, t])(r,0)=0.

We set some parameters and define the initial displacement:

c=1;R=2;
f[r_]=1-r/R;

Here are some alternate initial displacements:

f[r_]=(1-r/R)^4 (4 r/R+1);
(*f[r_]=Cos[r π/(2 R)];*)
(*f[r_]=Piecewise[{{1,r<1/2},{0,r<R}}];*)

Plot[f[r],{r,0,R}]

RevolutionPlot3D[f[r],{r,0,R}]

We showed that, for general initial position f and general initial velocity g, the solution is of the form
	u(r,t) = (∑_(n=1)^∞)((Subscript[a, n])cos((Sqrt[(Subscript[λ, n])])c t)+(Subscript[b, n])sin((Sqrt[(Subscript[λ, n])])c t))(Subscript[J, 0])((Sqrt[(Subscript[λ, n])])r)
with
	(Subscript[a, n])=(((∫_0^R)f(r)(Subscript[J, 0])((Sqrt[(Subscript[λ, n])])r)rdr)/((∫_0^R)((((Subscript[J, 0])((Sqrt[(Subscript[λ, n])])r)))^2)rdr))  and  (Subscript[b, n])=(((∫_0^R)g(r)(Subscript[J, 0])((Sqrt[(Subscript[λ, n])])r)rdr)/((∫_0^R)((((Subscript[J, 0])((Sqrt[(Subscript[λ, n])])r)))^2)rdr)).
Since we assumed the initial velocity to be zero we have (Subscript[b, n])=0.

First, we know that the eigenvalues are given by
	(Subscript[λ, n])= (((((Subscript[z, n])/R)))^2),
where (Subscript[z, n]) is the n-th zero of the Bessel function (Subscript[J, 0]).

The Bessel function (Subscript[J, 0]):

Plot[BesselJ[0,r],{r,0,50}]

The zeros of the Bessel function look almost equally spaced, but they are not. However, their spacing approaches π. Here are the first 16 zeros from the graph above along wih their spacing:

Grid[{{Grid[Table[{i,BesselJZero[0,i]//N,BesselJZero[0,i]//N-0},{i,1}]]},{Grid[Table[{i,BesselJZero[0,i]//N,BesselJZero[0,i]-BesselJZero[0,i-1]//N},{i,2,16}]]}},Alignment->"."]

Now, back to our calculations:

The eigenvalues are given by

Lambda[n_]=N[(BesselJZero[0,n]/R)^2];
Table[Lambda[n],{n,1,16}]

We now compute the Fourier coefficients (Subscript[a, n]):

a[n_] = Integrate[f[r]BesselJ[0,Sqrt[Lambda[n]]r] r,{r,0,R}]/Integrate[BesselJ[0,Sqrt[Lambda[n]]r]^2 r,{r,0,R}];
Table[a[n],{n,1,10}]

In our setting, the N-th partial sum of the Fourier series solution of the wave equation is
	u(r,t) = (∑_(n=1)^N)(Subscript[a, n])cos((Sqrt[(Subscript[λ, n])])c t)(Subscript[J, 0])((Sqrt[(Subscript[λ, n])])r).

u[r_,t_,N_]=Sum[a[n] Cos[Sqrt[Lambda[n]]c t]BesselJ[0,Sqrt[Lambda[n]]r],{n,1,N}];
uplot=u[r,t,20];

We plot the (partial sum approximation to the) solution at time t=0.

RevolutionPlot3D[uplot/.t->0,{r,0,R}]

Some more plots at different times t:

Plot[uplot/.t->0.5,{r,0,R}]
Plot[uplot/.t->1.0,{r,0,R}]
Plot[uplot/.t->1.5,{r,0,R}]

And now a little movie:

Manipulate[RevolutionPlot3D[uplot/.t->tplot,{r,0,R},PlotRange->{All,All,{-1,1}}],{tplot,0,10}]

Normal modes and overtones:

The modes are given by the individual terms of the Fourier series solution u, i.e.,
	(Subscript[a, n])cos((Sqrt[(Subscript[λ, n])])c t)(Subscript[J, 0])((Sqrt[(Subscript[λ, n])])r).

mode[r_,t_,n_]=a[n] Cos[Sqrt[Lambda[n]]c t]BesselJ[0,Lambda[n]r];
mode1=mode[r,t,1];
mode2=mode[r,t,2];
mode3=mode[r,t,3];

We can animate the first few modes (note that the amplitude decreases and the frequency increases from one mode to the next):

Manipulate[RevolutionPlot3D[mode1/.t->tplot,{r,0,R},PlotRange->{All,All,{-1,1}}],{tplot,0,10}]

Manipulate[RevolutionPlot3D[mode2/.t->tplot,{r,0,R},PlotRange->{All,All,{-0.1,0.1}}],{tplot,0,10}]

Manipulate[RevolutionPlot3D[mode3/.t->tplot,{r,0,R},PlotRange->{All,All,{-0.05,0.05}}],{tplot,0,10}]

Solution of wave equation on a rectangle:

We consider the 2D wave equation with c=1:
	(Subscript[u, tt])(x,y,t) = (Subscript[u, xx])(x,y,t) +(Subscript[u, yy])(x,y,t),   0≤x≤1,  0≤y≤2,
	u(0,y,t)=u(1,y,t)=u(x,0,t)=u(x,2,t)=0,
	u(x,y,0)=f(x,y),   (Subscript[u, t])(x,y,0)=g(x,y).

We define the initial displacement:

f[x_,y_]=Piecewise[{{x,x≤1/2},{1-x,x≤1}}]*Piecewise[{{y,y≤1},{2-y,y≤2}}];

Plot3D[f[x,y],{x,0,1},{y,0,2}]

We can compute the Fourier coefficients

a[n_,m_]=2 Integrate[Integrate[f[x,y] Sin[m π y/2],{y,0,2}]Sin[n π x],{x,0,1}];
Grid[Table[a[n,m],{n,1,10},{m,1,10}]]

Here are the eigenvalues

Lambda[n_,m_]=(n π)^2 + (m π/2)^2;

The solution, truncated at N and M, respectively, is given by

u[x_,y_,t_,N_,M_]:=Sum[Sum[a[n,m] Cos[Sqrt[Lambda[n,m]]t]Sin[n π x]Sin[m π y/2],{n,1,N}],{m,1,M}]
uplot=u[x,y,t,10,20];

And here is an animation of the approximate solution using N=M=10:

Manipulate[Plot3D[uplot/.t->tplot,{x,0,1},{y,0,2},PlotRange->{All,All,{-1/2,1/2}}],{tplot,0,5}]

And finally some modes:

mode[x_,y_,t_,n_,m_]=a[n,m] Cos[Sqrt[Lambda[n,m]]t]Sin[n π x]Sin[m π y/2];
mode11=mode[x,y,t,1,1];
mode13=mode[x,y,t,1,3];
mode33=mode[x,y,t,3,3];
Plot3D[mode11/.t->0,{x,0,1},{y,0,2}]
Plot3D[mode13/.t->0,{x,0,1},{y,0,2}]
Plot3D[mode33/.t->0,{x,0,1},{y,0,2}]

Solution of wave equation on a disk (with zero initial displacement):

We consider the general wave equation on a disk of radius R
	(Subscript[u, tt]) = (c^2) (((Subscript[(r (Subscript[u, r])), r])/r) +((Subscript[u, θθ])/(r^2)))
subject to the boundary condition
	u(R,θ,t)=0
and initial conditions
	u(r,θ,0)=f(r,θ),
	(Subscript[u, t])(r,θ,0)=g(r,θ).
If the problem is circularly symmetric then the PDE simplifies to
	(Subscript[u, tt]) = (c^2) ((Subscript[(r (Subscript[u, r])), r])/r).
If, in addition, we assume that the initial displacement is zero, then the boundary and initial conditions become
	u(R,t)=0,
	u(r,0)=0,
	(Subscript[u, t])(r,0)=g(r).

We set some parameters and define the initial displacement:

c=1;R=2;
g[r_]=-Exp[-10 r^2];
(*g[r_]=-Piecewise[{{1,r<1/2},{0,r<R}}];*)

Plot[g[r],{r,0,R}]

RevolutionPlot3D[g[r],{r,0,R},PlotRange->All]

We showed that, for general initial position f and general initial velocity g, the solution is of the form
	u(r,t) = (∑_(n=1)^∞)((Subscript[a, n])cos((Sqrt[(Subscript[λ, n])])c t)+(Subscript[b, n])sin((Sqrt[(Subscript[λ, n])])c t))(Subscript[J, 0])((Sqrt[(Subscript[λ, n])])r)
with
	(Subscript[a, n])=(((∫_0^R)f(r)(Subscript[J, 0])((Sqrt[(Subscript[λ, n])])r)rdr)/((∫_0^R)((((Subscript[J, 0])((Sqrt[(Subscript[λ, n])])r)))^2)rdr))  and  (Subscript[b, n])=(((∫_0^R)g(r)(Subscript[J, 0])((Sqrt[(Subscript[λ, n])])r)rdr)/((∫_0^R)((((Subscript[J, 0])((Sqrt[(Subscript[λ, n])])r)))^2)rdr)).
Since we assumed the initial velocity to be zero we have (Subscript[b, n])=0.

First, we know that the eigenvalues are given by
	(Subscript[λ, n])= (((((Subscript[z, n])/R)))^2),
where (Subscript[z, n]) is the n-th zero of the Bessel function (Subscript[J, 0]).

The eigenvalues are given by

Lambda[n_]=N[(BesselJZero[0,n]/R)^2];
Table[Lambda[n],{n,1,16}]

We now compute the Fourier coefficients (Subscript[b, n]):

b[n_] = Integrate[g[r]BesselJ[0,Sqrt[Lambda[n]]r] r,{r,0,R}]/Integrate[BesselJ[0,Sqrt[Lambda[n]]r]^2 r,{r,0,R}];
Table[N[b[n]],{n,1,10}]

In our setting, the N-th partial sum of the Fourier series solution of the wave equation is
	u(r,t) = (∑_(n=1)^N)(Subscript[b, n])sin((Sqrt[(Subscript[λ, n])])c t)(Subscript[J, 0])((Sqrt[(Subscript[λ, n])])r).

u[r_,t_,N_]=Sum[N[b[n]] Sin[Sqrt[Lambda[n]]c t]BesselJ[0,Sqrt[Lambda[n]]r],{n,1,N}];
uplot=u[r,t,20];

We plot the (partial sum approximation to the) solution at time t=0.

RevolutionPlot3D[uplot/.t->0,{r,0,R}]

Some more plots at different times t:

Plot[uplot/.t->0.5,{r,0,R}]
Plot[uplot/.t->1.0,{r,0,R}]
Plot[uplot/.t->1.5,{r,0,R}]

And now a little movie:

Manipulate[RevolutionPlot3D[uplot/.t->tplot,{r,0,R},PlotRange->{All,All,{-1,1}}],{tplot,0,10}]