Special Functions

We have been using Bessel functions as though they are as familiar as Sine and Cosine!  Folks have questions about what they are and where they come from.  We are going to try to answer that along with how many similar named functions are there.   Well Mathematica is pretty encyclopaedic about such things.

?BesselJ

If this is not enough stuff NIST has a wondrous web site. 
	http://dlmf.nist.gov/
with lots more information.

Series[ BesselJ[n,x],{x,0, 5}]
Series[ BesselY[n,x],{x,0, 5}]

Bessel Function Computation:

Math/matlab/Mathematica

The NIST site is terse about how to compute Bessel functions!  It says recurrsion is a good deal.  For more details take a look at the very readable (because it is pretty old) article 
http://www.ams.org/journals/mcom/1959-13-066/S0025-5718-1959-0105794-5/S0025-5718-1959-0105794-5.pdf

It is pretty explicit.  Use the recurrsion 
	(Subscript[J, m-1])[x]+(Subscript[J, m+1])[x]  2 m/x (Subscript[J, m])[x]
to solve backwards from a large index M and with starting conditions
	(Subscript[J, M+1])[x]=0  and (Subscript[J, M])[x]=a
down to your target m and then normalize.

The Mathematica internal implementation notes are simple but far from detailed.  https://reference.wolfram.com/language/tutorial/SomeNotesOnInternalImplementation.html 
Search for Bessel!

According to 
http://stackoverflow.com/questions/7810086/calculate-bessel-function-in-matlab-using-jm1-2mjm-jm-1-formula

Matlab in BESSELJ uses a MEX interface to a Fortran library by D. E. Amos. 
Refs
D. E. Amos, “A subroutine package for Bessel functions of a complex argument and nonnegative order”, Sandia National Laboratory Report, SAND85-1018, May, 1985.

D. E. Amos, “A portable package for Bessel functions of a complex argument and nonnegative order”, Trans. Math. Software, 1986.

Mathworld 
http://mathworld.wolfram.com/BesselFunction.html
is interesting because it defines the functions using the recursion and then says “also” satisfies the ODE.

Taylor Series.

So what happens if you try to use a Taylor Series.

Clear[n,x,TS]
n=0;
order=12;
TS[x_]=Normal[Series[ BesselJ[n,x],{x,0, order}]]
Plot[
{ BesselJ[n,x],TS[x]},
{x,0, 8}]

LogPlot[
Abs[BesselJ[n,x]-TS[x]],
{x,0, 8}]
LogPlot[
((Abs[BesselJ[n,x]-TS[x]])/(Abs[BesselJ[n,x]])),
{x,0, 8}]

Timing

Lets just see how long Mathematica takes to compute various Bessel functions

n=12;
x=RandomReal[{0,123},(10^3)];
AbsoluteTiming[ 
BesselJ[n, x];
]
AbsoluteTiming[ 
BesselY[n, x];
]

n=123;
x=RandomReal[{0,123},(10^3)];
AbsoluteTiming[ 
BesselJ[n, x];
]

n=123+I;
x=RandomReal[{0,123},(10^3)];
AbsoluteTiming[ 
BesselJ[n, x];
]

Of course you probably have no idea if this is fast or slow.

x=RandomReal[{0,123},(10^3)];
AbsoluteTiming[ 
Cos[ x];
]

ODE

The ODE is 
	(x^2)y'' + x y' + ((x^2)- (n^2))y=0
for (Subscript[J, 0])[0]=1  while (Subscript[J, n])[0]=0 for any other real value of n.  Note: x=0 is a special singular point for the ODE.

Series Expansion

How do I see the series expansion for J.

Clear[x,y,n,r,ODE,a, SeriesSub]
order=4;
SeriesSub= {
y[x]->Sum[(Subscript[a, i]) (x^(i+r)), {i,0, order}],
y'[x]->Sum[(Subscript[a, i]) (i+r)(x^(i+r-1)), {i,0, order}],
y''[x]->Sum[(Subscript[a, i]) (i+r)(i+r-1)(x^(i+r-2)), {i,0, order}]
};
ODE=(x^2)y''[x] + x y'[x] + ((x^2)- (n^2))y[x];
TableForm[
CoefficientList[
(1/(x^r))ODE/.SeriesSub,
x],
TableHeadings->Automatic
]

Orthogonality

n=0;
{k1,k2}={4,6};
{α1,α2}=N[BesselJZero[n,{k1,k2}]]
Plot[ {
BesselJ[n,α1 x],
BesselJ[n,α2 x]
},{x,0, 1}]

Plot[ {
x BesselJ[n,α1 x]*BesselJ[n,α2 x]
},{x,0, 1}]
NIntegrate[ 
x BesselJ[n,α1 x]*BesselJ[n,α2 x],{x,0, 1}]

Legendre Computation:

Series Definition from ODE

Clear[x,y,n,r,ODE,a, SeriesSub]
order=4;
SeriesSub= {
y[x]->Sum[(Subscript[a, i]) (x^(i+r)), {i,0, order}],
y'[x]->Sum[(Subscript[a, i]) (i+r)(x^(i+r-1)), {i,0, order}],
y''[x]->Sum[(Subscript[a, i]) (i+r)(i+r-1)(x^(i+r-2)), {i,0, order}]
};
ODE=(1-(x^2)) y''[x]- 2 x y'[x]+n(n+1) y[x]
TableForm[
CoefficientList[(ODE/(x^(r-2)))/.SeriesSub,x],
TableHeadings->Automatic
]