polarVibes_BesselFns

This notebook contains code to illustrate Bessel functions and their application to solving problems involving the Laplacian inpolar coordinates.

Once again, define a style for plots so they project visibly.

fat = AbsoluteThickness[3]

Bessel functions

What does Bessel's equation (which often arises from separating the radial dependence of the Laplacian) look like and what are its solutions?

DSolve[x^2*y''[x]+x*y'[x]+(x^2-n^2)y[x]==0, y[x], x]

For each n, there is a part proporational to BesselJ[n,x], a Bessel function of the first kind of order n, and a part proportional to BesselY[n,x], a Bessel function of the second kind of order n.
What do these look like?

Table[Plot[{BesselJ[n,x], BesselY[n,x]}, {x,0,25}, PlotStyle -> fat, ImageSize -> 400, AxesLabel -> {x, {Style[(Subscript[J, n])[x], FontColor-> Blue, FontSize -> 24], Style[(Subscript[Y, n])[x], FontColor-> Purple, FontSize -> 24]}}], {n,0,3}]

So in each pair,  BesselY[n,x]  blows up at the origin (and so does not contribute to the solution on a full disk).

With each unit increment of n, BesselJ[n,x] has one more derivative that vanishes at the origin.

In each case, the functions asymptotically oscillate (so infinitely many roots are guaranteed) and decay like (Sqrt[x]).

Table[Plot[{x*BesselJ[n,x]^2, x*BesselY[n,x]^2}, {x,0,25}, PlotStyle -> fat, ImageSize -> 400, AxesLabel -> {x, {Style[(Subscript[J, n])[x]^2, FontColor-> Blue, FontSize -> 24], Style[(Subscript[Y, n])[x]^2, FontColor-> Purple, FontSize -> 24]}}], {n,0,3}]

Use Bessel functions as eigenfunctions for solving PDEs.

In the example of the circular drumhead, the ODE for the radial dependence turned out to be Bessel's equation where the order n is determined by the periodicity requirement in the θ direction. Specifically, if the trig functions in the θ direction have n cycles around the drum (i.e. in 2π radians), then the relevant radial solutions are (Subscript[J, n])(kr). (Remember that the other solution, (Subscript[Y, n])(kr)would blow up at r=0, so it cannot contribute to describe the finite displacement of the drumhead.)

For the moment, let's simplify things a bit and consider initial conditions that are independent of θ, so the solution becomes independent of θ. With this assumption, n=0 and the angular solution solution simplifies to Θ(θ)~ Cos[0*θ]~1. The time dependence is harmonic with frequency c*k, but we still do not know what k is. We need to determine it from the radial boundary condition, i.e. R(ρ) = 0, which implies (Subscript[J, 0])(kρ)=0. The eigenfunctions are then Bessel functions with arguments chosen so that the outer rim of the drumhead corresponds to a zero of the Bessel function.To have a hope of fitting arbitrary functions, we need to have an infinite number of eigenfunctions and an infinite number of associated roots of the Bessel function. The good news is that the Bessel functions oscillate (it can be shown that for large values of the argument, (Subscript[J, n])(x)~(Sqrt[(2/(πx))])cos(x-(π/4)-((nπ)/2))) and provide an infinite number of zeros. Are they evenly spaced like the zeros of sin and cos? Not quite (although they become evenly spaced for large arguments), so they will need to be found numerically. The good news is that this has been done for us automatically by BesselJZero. Use it to construct a list of the first several eigenvalues:

α = Table[BesselJZero[0,m],{m,1,5}]  //N

%-RotateRight[%]

So now we have the first 5 roots of (Subscript[J, 0])(x).Let's pick a value, say ρ=1, for the radius of the drumhead so we can generate some plots numerically. Radial solutions then look like (Subscript[J, 0])((Subscript[α, n])r). Let's take a look at these functions.

Table[BesselJ[0, α[[n]]*r], {n,1,5}]

Plot[Evaluate[%], {r,0,1},PlotStyle -> fat]

So we now have a sequence of radial functions that satisfy the homogeneous boundary condition, but do they provide an orthogonal basis? Let's look at the inner products (rembering the weighting factor r in the inner product that we obtain from Sturm-Liouville):

Table[
Integrate[
r*BesselJ[0, α[[m+1]]*r]*BesselJ[0, α[[n+1]]*r],
{r,0,1}],
{m,0,4}, {n,0,4}];

MatrixForm[Chop[%]]

So we have verified the orthogonality relation and we are in position to deal with finding coefficients to satisfy initial conditions. Again let's start simple and let the initial displacement be u(r,0)=(1-r) with zero initial velocity. Our approximate solution (including the first 5 terms) looks like:

c = 1;

u[r_, t_] := Sum[a[n]*Cos[c*α[[n]]*t]*BesselJ[0, α[[n]]*r], {n,1,5}]

u[r,0]

How do we determine the a[n] coefficients? Project the initial displacement onto the orthogonal set of Bessel functions:

a[n_] := N[Integrate[(1-r)*BesselJ[0,α[[n]]*r]*r,
{r,0,1}]/Integrate[BesselJ[0,α[[n]]*r]^2*r,
{r,0,1}]]

Table[a[n], {n,1,5}]

u[r,0]

Plot[
Evaluate[u[r,0]], {r,0,1}, PlotRange -> {0,1}, PlotStyle -> fat]

Again, this does not perfectly reproduce the specified initial displacement, but it is pretty good with just 5 terms included.

What do the vibrations look like? Let's evaluate at a sequence of times and plot each result. The collection of graphics can then be animated (or "GraphicsGrid"-ed).

cone = Table[
ParametricPlot3D[
Evaluate[{r*Cos[θ],r*Sin[θ],u[r,t]}],
{r,0,1}, {θ,0,2Pi},  PlotRange -> {-1,1}],
{t,0,2.3,.1}];

ListAnimate[cone]

Show[GraphicsGrid[ Partition[cone,3] ] , ImageSize -> 400]

Alternatively, we can think of this solution as a linear combination of modes (corresponding to the terms in our solution). What do the modes look like?

First Mode

Cos[c*α[[1]]*t]*BesselJ[0, α[[1]]*r]

mode1 = Table[
ParametricPlot3D[
Evaluate[{r*Cos[θ],r*Sin[θ],Cos[c*α[[1]]*t]*BesselJ[0, α[[1]]*r]}],
{r,0,1}, {θ,0,2Pi},  PlotRange -> {-1,1}],
{t,0,2.3,.1}];

ListAnimate[mode1]

Second mode

mode2 = Table[
ParametricPlot3D[
Evaluate[{r*Cos[θ],r*Sin[θ],Cos[c*α[[2]]*t]*BesselJ[0, α[[2]]*r]}],
{r,0,1}, {θ,0,2Pi},  PlotRange -> {-1,1}],
{t,0,2.3,.1}];

ListAnimate[mode2]

Third Mode

mode3 = Table[
ParametricPlot3D[
Evaluate[{r*Cos[θ],r*Sin[θ],Cos[c*α[[3]]*t]*BesselJ[0, α[[3]]*r]}],
{r,0,1}, {θ,0,2Pi},  PlotRange -> {-1,1}],
{t,0,2.3,.1}];

ListAnimate[mode3]

So each higher mode has an additional nodal line (where the displacement is zero) and a higher frequency.

Include ϑ dependence

What happens if you do not assume axial symmetry?

Let's look at the modes. Let n be the frequency of the ϑ dependence and let m denote which root of the (n^th) Bessel function we are using. How should we order the modes?

zero = Table[BesselJZero[n,m],{n,0,5},{m,1,5}] //N

Note that the n+1 below accounts for the fact that the index of an array starts at 1 (instead of zero, the initial index for the Bessel functions).

mode[n_, m_, r_, θ_, t_] := Cos[c*zero[[n+1,m]]*t]*
Cos[n*θ]* 
BesselJ[n, zero[[n+1,m]]*r]

mode[0,1,r,θ,0] /. {r->0.5,θ->Pi}

Plot[mode[0,1,r,θ,0] /. θ->0, {r,0,1}]

ParametricPlot3D[Evaluate[{r*Cos[θ], r*Sin[θ], mode[0,1,r,θ,0]}], {r,0,1}, {θ,0,2Pi}]

ParametricPlot3D[{r*Cos[θ], r*Sin[θ], mode[1,1,r,θ,0]}, {r,0,1}, {θ,0,2Pi}]

ParametricPlot3D[{r*Cos[θ], r*Sin[θ], mode[2,1,r,θ,0]}, {r,0,1}, {θ,0,2Pi}]

ParametricPlot3D[{r*Cos[θ], r*Sin[θ], mode[0,2,r,θ,0]}, {r,0,1}, {θ,0,2Pi}]

ParametricPlot3D[{r*Cos[θ], r*Sin[θ], mode[1,2,r,θ,0]}, {r,0,1}, {θ,0,2Pi}]

These start to get a bit tough to decipher, so create animations (as above) or create ContourPlots (including nodal patterns).

period22 = 2Pi/zero[[2,1]]

frames22 = Table[ParametricPlot3D[{r*Cos[θ], r*Sin[θ], mode[1,2,r,θ,t]}, {r,0,1}, {θ,0,2Pi},
PlotRange -> {{-1,1},{-1,1},{-0.6,0.6}}],
{t,0,period22, period22/15}];

ListAnimate[frames22]

Semi-annular drum

Characteristic equation: (Subscript[J, n])[βa](Subscript[Y, n])[2βa]=(Subscript[J, n])[2βa](Subscript[Y, n])[βa]

Plot[Evaluate[{BesselJ[n,2x]*BesselY[n,x],
BesselJ[n,x]*BesselY[n,2x]} /. n->0], {x,0,25}, PlotStyle -> fat]

This looks like a mess, so let's plot the difference and look for behavior of the roots.

char[x_, n_] := BesselJ[n,2x]*BesselY[n,x]-
BesselJ[n,x]*BesselY[n,2x]

Plot[Evaluate[Table[char[x,n],{n,1,2}]], {x,0,10}]

Now we see a nice oscillation (suggesting and infinite set of roots) with  asymptotically regular spacing.

Consider a particular mode: n=1, m=1

β[1,1] = FindRoot[char[x,1], {x,3}][[1,2]]

f[r_, 1,1] :=
BesselJ[ 1,β[1,1]*r]- (BesselJ[1,β[1,1]]/BesselY[1,β[1,1]])*BesselY[ 1,β[1,1]*r]

Plot[f[r,1,1],{r,1,2}]

ParametricPlot3D[{r*Cos[θ], r*Sin[θ],f[r,1,1]*Sin[θ]},{r,1,2}, {θ,0,Pi}, AxesLabel -> {X,Y,U}, BoxRatios ->{2,1,0.5}]

Consider a particular mode: n=2, m=1

β[2,1] = FindRoot[char[x,2], {x,3}][[1,2]]

f[r_, 2,1] :=
BesselJ[ 2,β[2,1]*r]- (BesselJ[2,β[2,1]]/BesselY[2,β[2,1]])*BesselY[ 2,β[2,1]*r]

Plot[f[r,2,1],{r,1,2}]

ParametricPlot3D[{r*Cos[θ], r*Sin[θ],f[r,2,1]*Sin[2θ]},{r,1,2}, {θ,0,Pi}, AxesLabel -> {X,Y,U}, BoxRatios ->{2,1,0.5}]

Consider a particular mode: n=1, m=2

β[1,2] = FindRoot[char[x,1], {x,3+Pi}][[1,2]]

f[r_, 1,2] :=
BesselJ[ 1,β[1,2]*r]- (BesselJ[1,β[1,2]]/BesselY[1,β[1,2]])*BesselY[ 1,β[1,2]*r]

Plot[f[r,1,2],{r,1,2}]

ParametricPlot3D[{r*Cos[θ], r*Sin[θ],f[r,1,2]*Sin[θ]},{r,1,2}, {θ,0,Pi}, AxesLabel -> {X,Y,U}, BoxRatios ->{2,1,0.5}]

Check orthogonality...