Bessel Functions
BesselJ[n,x]

by
K. D.  Stroyan
University of Iowa

copyright 1997 by Academic Press, Inc. - All rights reserved.

Special Mathematica Functions used in this NoteBook

Off[General::spell1];
Needs["Graphics`Colors`"]
seriesGraphs[f_, series_, radius_, numterms_Integer, opns___Rule] :=
	Module[ {functions, graylevels, i, j,n,x,terms},
	Do[
	functions = Table[series[x,n], {n, 0, i}] ;
	AppendTo[functions, f[x]] ;
	graylevels = Table[{GrayLevel[j]}, {j, 6/10, 0, -6/(10 i)}] ; 
	graylevels=AppendTo[graylevels,{AbsoluteThickness[2],Red}];
	Plot[Release[functions], {x, -radius, radius},
		PlotStyle -> graylevels,PlotLabel ->i"terms",
		opns,AspectRatio->Automatic];
	,{i, 1, numterms}];
	Do[
	Plot[Release[{f[x], series[x,n]}],{x, -radius, radius},
		PlotLabel ->n"terms ", opns,
		PlotStyle -> {{AbsoluteThickness[2],Red},{GrayLevel[0]}},
		AspectRatio->Automatic];
	,{n, 1, numterms}];
	];

Needs["Graphics`Colors`"]
powerGraphs[f_,  radius_, numterms_, opns___Rule] :=
	Module[ {functions, graylevels, i, j,n,x,terms},
	Do[
	functions = Table[Normal[Series[f[x],{x,0,n}]], {n, 0, i}] ;
	AppendTo[functions, f[x]] ;
	graylevels = Table[{GrayLevel[j]}, {j, 6/10, 0, -6/(10 i)}] ; 
	graylevels=AppendTo[graylevels,{AbsoluteThickness[2],Red}];
	Plot[Release[functions], {x, -radius, radius},
		PlotStyle -> graylevels, PlotLabel ->i"terms",
		opns,AspectRatio->Automatic];
	,{i, 1, numterms}];
	Do[
	Plot[Release[{Normal[Series[f[x],{x,0,n}]],f[x]}],{x, -radius, radius},
		PlotLabel ->n"terms ", opns, 
		PlotStyle -> {{GrayLevel[0]},{AbsoluteThickness[2],Red}},AspectRatio->Automatic];
	,{n, 1, numterms}];
	];
On[General::spell1];

Notebook Overview: Exercise 26.4.3

Some functions that arise in mathematics or science are known mainly by their Taylor series or by a differential equation.  Famous examples are the Bessel functions which are similar to the Sin and Cos functions in some respects and arise in many applications through their differential equation.  This differential equation comes from the heat, wave or Laplace partial differential equations in radially symmetric problems.  Substituting polar or cylindrical coordinates leads to an ordinary differential equation that, by definition, describes the Bessel function.

Bessel's equation (with parameter n) is:

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

One solution of this equation with is called "J[n,x]" or BesselJ[n,x] in built-in Mathematicaese.  Mathematica  can work directly with BesselJ[n,x].

\(n\  = \ 0; 

    ParametricPlot3D[{r\ Cos[θ], \ r\ Sin[θ], \ 
        BesselJ[n, r]}, {r, 0, 15}, {θ, \(-π\), π}, 
      PlotPoints -> 30, BoxRatios -> {2, 2, 1}]; \)

\(n\  = \ 1; 

    ParametricPlot3D[{r\ Cos[θ], \ r\ Sin[θ], \ 
        BesselJ[n, r]}, {r, 0, 15}, {θ, \(-π\), π}, 
      PlotPoints -> 30, BoxRatios -> {2, 2, 1}]; \)

References to the Text

This NoteBook goes with Exercise 26.4.3 of the text.

A First Look

Mathematica can find power series expansions of functions.  <Enter> the next example.

Mathematica's Built-in Series[.] function

<Enter> the next computation to see Mathematica's description of the Series[.] function.

?? Series

BesselJ[n, x]

m = 10; Print[Series[BesselJ[0, x], {x, 0, m}]]; Print[Series[BesselJ[1, x], {x, 0, m}]];

Comparison of Formulas

The formulas for these terms given in Exercise 26.4.3 is as follows:

Clear[j0, j1, x]; j0[x_] = 
      (∑_(k = 0)^5)
          (((((-1))^k) (x^(2 k)))/(((k!)^2) (2^(2 k))))j1[x_] = 
      (∑_(k = 0)^5)
          (((((-1))^k) (x^(2 k + 1)))/(k! (k + 1)! 
            (2^(2 k + 1))))

Compare these results with the results of the Series[.] command above.

The Series Graphs

The powerGraphs[.] function defined above in the Initialization cell uses the built-in Series[.] function to plot a function and its power series approximations within the "radius" you specify and up to the number of terms "n" you specify.  (The graphs only converge inside the "radius of convergence" studied in Chapter 26.)

\(radius = 15; 
f[x_]\  = \ BesselJ[0, x]; 

    powerGraphs[f, radius, m, AspectRatio -> 1, 
      PlotRange -> {\(-\(radius\/5\)\), radius\/5}]; \)

\(radius = 15; 
f[x_]\  = \ BesselJ[1, x]; 

    powerGraphs[f, radius, m, AspectRatio -> 1, 
      PlotRange -> {\(-\(radius\/5\)\), radius\/5}]; \)

Practice Exercise

Use Mathematica to find the power series expansion of BesselJ[2,x] and plot the first 20 terms.

Exercise 26.4.3: The Differential Equations and Series

(0) Uniqueness?

Solutions to initial value problems are unique provided the growth law is smooth.  Look up the exact statement of the Existence & Uniqueness Theorem from Chapter 22.

Convert the Besssel Equation

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

into a first order system of the form:

y'[x] == f[y[x],z[x]]
	z'[x] ==  g[y[x],z[x]]

And say which initial conditions make it satisfy the theorem. In particular, show that the theorem does NOT apply to:

y[0] == y0,
	z[0] == z0,
	y'[x] == f[y[x],z[x]],
	z'[x] ==  g[y[x],z[x]]

(HINT: What are the functions  f[y,z] and g[y,z]  in the formulas?  Where are they smooth?)

(1) BesselJ[n,x] Series Solutions

Simple Verification for BesselJ[0,x]

Using the textbook formual for the series of BesselJ[0,x], we can verify solution of the IVP as follows:

The formula up to (x^10)

Clear[m, j0]; m = 5; j0[x_] = 
      (∑_(k = 0)^m)
          (((((-1))^k) (x^(2 k)))/(((k!)^2) (2^(2 k))))

The initial values

j0[0]D[j0[x], x] /. x -> 0

Substitution in the differential equation

\(x\^2\ D[j0[x], {x, 2}]\  + \ x\ D[j0[x], x]\  + \ x\^2\ j0[x]
\)Simplify[ 
      (x^2) D[j0[x], {x, 2}]  +  x D[j0[x], x]  +  (x^2) j0[x]]

Why is the (x^12)  term left over?  (What would happen if you defined j0[x] as the sum to m = 20?)

Series solution from scratch:

Solve
	y[0] == 1,
	y'[0] == 0,
	(x^2) y''[x] + x y'[x] +((x^2)) y[x] == 0
	
	by substituting the unknown series
	
	y[x] = a[0] + a[1] x + a[2] (x^2) + ...+ a[n] (x^n) + ...
	
	y'[x] = a[1] + 2 a[2] x + 3 a[3] (x^2) + ...+ (n+1) a[n+1] (x^n) + ...
	
	y''[x] = 2 a[2] + 3*2 a[3] x + ... + (n+2)(n+1) a[n+2] (x^n) + ...
	
	into the initial conditions and differential equation.  Here's a start:
	
	INITIAL CONDITIONS:
	y[0] = a[0] = 1
	y'[0] = a[1] = 0
	
	THE EQUATION:
		(x^2) y''[x] + x y'[x] + (x^2) y[x] == 0
		
2 a[2] (x^2) + 3*2 a[3] (x^3) + 4*3*a[4] (x^4) +  ... + (n)(n-1) a[n] (x^n) + ...
+
0 x + 2 a[2] (x^2) + 3 a[3] (x^3) + ...+ (n) a[n] (x^n) + ...
+
 1 (x^2) + 0 (x^3) + a[2] (x^4) + ...+ a[n-2] (x^n) + ...
=
0

So
	2a[2] + 2a[2] + 1 == 0 => a[2] = -1/4, 
	a[3] = 0,
	4*3 a[4] + 4 a[4] + a[2] == 0, => a[4] = 1/((2^2)*(2^4)),

aODD = a[2n+1] = 0
	
	solutionJ[o, x] = 
        1 - ((x^2)/4) + ((x^4)/64) + ⋯ + 
          (((((-1))^n) (x^(2 n)))/(((n!)^2) (2^(2 n)))) + 
          ⋯

Automatic series solution

Clear[n, a]; series = (∑_(k = 0)^10) a[k] (x^k) + O[x]^10

t0  =  (x^2  -  n^2)*seriest1  =  x*D[series,  x]t2  =  x^2*D[series,  {x,  2}]

Solve for Ten Coefficients of BesselJ[0,x]:

a[0] = 1; a[1] = 0; n = 0; Solve[t2 + t1 + t0 == 0]

Compare coefficients:

Series[BesselJ[0, x], {x, 0, 8}]

(2) BesselJ[1,x] Series Solution

Solve

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

by one of the three series approaches above for BesselJ[0,x], showing that the series:

j1[x] = 
        (x/2) - ((x^3)/4) + ((x^4)/64) + ⋯ + 
          (((((-1))^n) (x^(2 n)))/(((n!)^2) (2^(2 n)))) + 
          ⋯

satisfies the initial value problem.

(3) Convergence

Use the ratio test to show that the series

j0[x]

and

j1[x]

both converge uniformly for  | x | ≤ r, for any  r > 0. (The radius of convergence is +∞.)

(4) Identities via Series

Use the series to show that

j0'[x] = -j1[x]

and

x j0[x] = d(x j1[x])/dx

j0[x]j1[x]j0'[x]+j1[x]

(6) Differentiate Analytically

\(Clear[n]; 
D[BesselJ[n, x], x]\)

(7) Analytical Solutions

BesselJ[n,x]  is a solution to the initial value problem:

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

Humm, HyperGeometric functions...  Another new kind of function.

Let's try the phase variable trick:

Clear[x, y, z, n]; DSolve[{y'[x]==z[x],(x^2) z'[x]==((n^2) - (x^2)) y[x] - x z[x]},{y[x], z[x]},x]

Notice that Mathematica is confused and didn't find THE solution.

BesselY[n,x]: Another Solution to the Non-unique Solution

ONE NON-unique solution to this problem also goes by the name BesselJ[n,x]. It does satisfy the equation:

\(Clear[n]; \ 
BesselJ[n, \ 0]\)D[BesselJ[n,  x],  x]

D[BesselJ[n,  x],  {x,  2}]

\(x\^2*D[BesselJ[n, \ x], \ {x, \ 2}]\  + \ x*D[BesselJ[n, \ x], \ x]\  + 
      \ 
\ \ \((x\^2 - \ n\^2)\)*BesselJ[n, \ x]\)

\(n\  = \ 0; 

    x\^2*D[BesselJ[n, \ x], \ {x, \ 2}]\  + \ x*D[BesselJ[n, \ x], \ x]\  + \ 
      
\ \ \((x\^2 - \ n\^2)\)*BesselJ[n, \ x]\)

This is not zero, but if  n  is an integer, there are relations between the various Bessel functions.  For example,

\(n\  = \ 1; 

    Plot[{BesselJ[n, x], BesselJ[\(-n\), x]}, {x, 0, 20}, 
      PlotStyle -> {Red, Blue}]; 
Series[BesselJ[n, x], {x, 0, 10}]\)Series[BesselJ[-n, x], {x, 0, 10}]

\(n\  = \ 2; 

    Plot[{BesselJ[n, x], BesselJ[\(-n\), x]}, {x, 0, 20}, 
      PlotStyle -> {{Thickness[0.02], Red}, Blue}]; 

    Series[BesselJ[n, x], {x, 0, 10}]\)Series[BesselJ[-n, x], {x, 0, 10}]

Using these identities we see that the BesselJ functions satisfy the differential equation.  In terms of the series:

\(n = 0; 

    Series[x\^2*D[BesselJ[n, \ x], \ {x, \ 2}]\  + \ 
        x*D[BesselJ[n, \ x], \ x]\  + \ 
\ \ 
        \((x\^2 - \ n\^2)\)*BesselJ[n, \ x], {x, 0, 10}]\)

\(n = 1; 

    Series[x\^2*D[BesselJ[n, \ x], \ {x, \ 2}]\  + \ 
        x*D[BesselJ[n, \ x], \ x]\  + \ 
\ \ 
        \((x\^2 - \ n\^2)\)*BesselJ[n, \ x], {x, 0, 10}]\)

Information["BesselJ", LongForm -> True]

Information["BesselY", LongForm -> False]

Clear[y, x]; y[x_] := BesselY[1, x]; n = 5; Print[Series[y[x], {x, 0, n}]]; radius = 15; powerGraphs[y, radius, n, AspectRatio -> 1, 
      PlotRange -> {-(radius/5), (radius/5)}];

N[EulerGamma]