Bessel Functions

Bessel's Equation

Here is Bessel's equation:

L[y_]:=(x^2)D[y,{x,2}]+x*D[y,{x,1}]+((x^2)-(v^2))y==0
L[y[x]]

Mathematica can solve this, and we get:

DSolve[L[y[x]],y[x],x]

Our goal here is to understand what how these two solutions are obtained. The solutions are Bessel Functions of the First Kind of order v ((Subscript[J, v])(x)) and Bessel Functions of the Second Kind of order v ((Subscript[Y, v])(x)).

Bessel Functions of The First Kind of Order v, (Subscript[J, v])(x)

The handout explains in detail how the Bessel Function of the First Kind is arrived at as the series solution to Bessel's equation for the root of the indicial equation r=v>0. Contained here are simply some mathematical exploratiosn along the way, and plots of the Bessel function (Subscript[J, v])(x).

f[m_]:=(∏_(k=1)^m)(k+v)
f[m]

(Subscript[y1, v_])[x_]=((((x/2)))^v)(∑_(m=0)^∞)(((((-1))^m)(x^(2m)))/((2^(2m))m! Gamma[1+v+m]))

The Bessel function reduces to a sine function is certain cases.

(Subscript[y1, 1/2])[x]

The Bessel function is defined for complex values.

(Subscript[y1, -I])[2.0]

Here is a sketch of (Subscript[J, v])(x).

functionlist=Table[(Subscript[y1, v])[x],{v,0,5}];
ps1=Table[{Dashing[{0.02,k/5*0.03}],Thick},{k,0,5}];
Plot[Evaluate[functionlist],{x,0,10},PlotStyle->ps1,ImageSize->{500,Automatic},PlotRange->{Automatic,{-1,1}},AspectRatio->1/3,BaseStyle->20]

Bessel Functions of The Second Kind of Order v, (Subscript[Y, v])(x)

Finding the second solution to Bessel's equation is a complicated process in general. Here, we look at finding the first few terms in the second solution using reduction of order.  Here we use the formula for the second solution we found when we studied reduction of order in general earlier.

y2(x)=(∫((e^(-∫p(x)dx))/(((y1(x)))^2))dx) y1(x)

The most difficult aspect of the calculation is the squaring and inversion of a series, which is easy to do with a computer algebra system. This might not always be feasible! Here we are lucky in that we have the first solution as a Bessel function, so we do not have to work with the infinite series at all.

(Subscript[y2, v_])[x_]=(∫((E^(-∫(1/x)\[DifferentialD]x))/((((Subscript[y1, v])[x]))^2))\[DifferentialD]x) *(Subscript[y1, v])[x]

Notice that we get the second solution as a Bessel Function of the Second Kind of order v, up to a constant multiple. To do this integral by hand we would need to use a Big Book of Integrals, or spend more time learning about Bessel functions.

functionlist=Table[(2/π)(Subscript[y2, v])[x],{v,0,5}];
ps1=Table[{Dashing[{0.02,k/5*0.03}],Thick},{k,0,5}];
Plot[Evaluate[functionlist],{x,0,10},PlotStyle->ps1,ImageSize->{500,Automatic},PlotRange->{Automatic,{-1,1}},AspectRatio->1/3,BaseStyle->20]

Application: Vibrating Drumhead (Looking forward to PDEs)

Here is the wave equation in polar coordinates, which is a partial differential equation (PDE). So this is looking ahead to how the Bessel equation arises when solving certain partial differential equations. The details of how all this works would be in a PDE class, which you are now ready to take, having completed an ODE course!

For PDEs, we need more than just initial conditions, we usually have some boundary conditions as well. Here is the PDE and conditions for a vibrating circular drumhead in polar coordinates.

Here is the solution, without any explanation of where it comes from. Technique: Separation of Variables, which you would learn in a PDE class.

a=1;
c=1;
f[r_,θ_]=(1-(r^2))r Sin[3θ];
g[r_,θ_]=(1-(r^2))(r^2) Sin[2θ];

(Subscript[λ, m_,n_]):=((BesselJZero[m,n])/a)//N

(Subscript[a, 0,n_]):=(2/( π (a^2)((BesselJ[1,a (Subscript[λ, 0,n])])^2) ))(∫_0^a)(∫_0^(2π))f[r,θ]BesselJ[0,r (Subscript[λ, 0,n])]r\[DifferentialD]θ\[DifferentialD]r
(Subscript[a, m_,n_]):=(2/( π (a^2)((BesselJ[m+1,a (Subscript[λ, m,n])])^2) ))(∫_0^a)(∫_0^(2π))f[r,θ]Cos[m θ]BesselJ[m,r (Subscript[λ, m,n])]r\[DifferentialD]θ\[DifferentialD]r
(Subscript[b, m_,n_]):=(2/( π (a^2)((BesselJ[m+1,a (Subscript[λ, m,n])])^2) ))(∫_0^a)(∫_0^(2π))f[r,θ]Sin[m θ]BesselJ[m,r (Subscript[λ, m,n])]r\[DifferentialD]θ\[DifferentialD]r

(Subscript[a2, 0,n_]):=(2/(c (Subscript[λ, 0,n]) π (a^2)((BesselJ[1,a (Subscript[λ, 0,n])])^2) ))(∫_0^a)(∫_0^(2π))g[r,θ]BesselJ[0,r (Subscript[λ, 0,n])]r\[DifferentialD]θ\[DifferentialD]r
(Subscript[a2, m_,n_]):=(2/(c (Subscript[λ, m,n]) π (a^2)((BesselJ[m+1,a (Subscript[λ, m,n])])^2) ))(∫_0^a)(∫_0^(2π))g[r,θ]Cos[m θ]BesselJ[m,r (Subscript[λ, m,n])]r\[DifferentialD]θ\[DifferentialD]r
(Subscript[b2, m_,n_]):=(2/(c (Subscript[λ, m,n]) π (a^2)((BesselJ[m+1,a (Subscript[λ, m,n])])^2) ))(∫_0^a)(∫_0^(2π))g[r,θ]Sin[m θ]BesselJ[m,r (Subscript[λ, m,n])]r\[DifferentialD]θ\[DifferentialD]r

sol[r_,θ_,t_,k_]:=(∑_(m=0)^k)(∑_(n=1)^k)BesselJ[m,r (Subscript[λ, m,n])]((Subscript[a, m,n])Cos[m θ]+(Subscript[b, m,n])Sin[m θ])Cos[c (Subscript[λ, m,n]) t]+(∑_(m=0)^k)(∑_(n=1)^k)BesselJ[m,r (Subscript[λ, m,n])]((Subscript[a2, m,n])Cos[m θ]+(Subscript[b2, m,n])Sin[m θ])Sin[c (Subscript[λ, m,n]) t]

solution[r_,θ_,t_]=sol[r,θ,t,4]

Notice those Bessel functions! The separation of variables leads to a Bessel equation ODE which has to be solved.

Notice the solution relies on knowing the location of the zeros of the Bessel functions.

Notice the solution is represented as an infinite series of sines and cosines, with coefficients that are in terms of integrals over Bessel functions (the a’s and b’s). The infinite series of sines and cosines is a Fourier series, also studied in a PDE class.

Notice I have to truncate the series to get an approximation to the solution.

Here is a plot of the solution as a function of time.

step=1/20;
endt=8;
Do[

plot1[t]=RevolutionPlot3D[solution[r,θ,t],{r,0,1},{θ,0,2Pi},PlotRange->{{-1,1},{-1,1},{-1.0,1.0}},Boxed->False,Axes->False,ImageSize->{600,Automatic}],

{t,0,endt,step}]

Manipulate[Show[plot1[t]],{t,0,endt,step}]

If you slow this one down, it shows the motion of the drum head quite well, in terms of where the extrema are located and how they move around.

step=1/50;
endt=3;
Do[

plot2[t]=RevolutionPlot3D[solution[r,θ,t],{r,0,1},{θ,0,2Pi},PlotRange->{{-1,1},{-1,1},{-1.0,1.0}},Boxed->False,Axes->False,ImageSize->{600,Automatic},ColorFunction->Function[{x,y,z,t,θ,r},Hue[z]]],

{t,0,endt,step}]

Manipulate[Show[plot2[t]],{t,0,endt,step}]