Bessel function

Bessel function referes to a solution $f$ of the Bessel equation

$ \!\!\!\!\!\!\!\!\! (1) ~ ~ ~ f(z)+f'(z)/z+(1-\nu/z^2)f(x) =0$

with the following asymptotic behaviour:

$\displaystyle f(x) \approx x^\nu \left( \frac{2^{-\nu}}{\mathrm{Factorial}(\nu)}+ O(x^2) \right)$

In particular,

For $\nu=0$, the initial condition is assumed $~f(0)=1~$, $~f'(0)=0$ , and

For $\nu=1$, the initial condition is assumed $~f(0)=0~$, $~f'(0)=1$;

Such a solution $f$ denoted with $\mathrm{BesselJ}[\nu,x]$ or with $J_\nu(x)$.

BesselJ

Due to singularity of the equation at $z=0$, the regular solution should have specific behavior. This solution is called $J_\nu$. For $\nu=0$ and $\nu=1$, there are specific implementations BesselJ0 and BesselJ1. Many formulas about the Bessel functions below are borrowed from the handbook by Abramowirtz,Stegun ; the numeration of formulas from there is used below.

Integral representations

$ \!\!\!\!\!\!\!\!\!\! (9.1.20) ~ ~ ~ \displaystyle

J_\nu(z) = \frac{(z/2)^{\nu}}{\pi^{1/2} ~(\nu-1/2)!} ~ \int_0^\pi ~ \cos(z \cos(t)) \sin(t)^{2 \nu} ~t~ \mathrm d t $

Sonin representation

The Mehler,Sonin formulas suggest that

$\displaystyle J_\nu(z)=\frac{2}{\pi} \int_0 ^\infty \sin(x \cos(t) - \pi \nu/2) \cos(\nu t) \mathrm d t$

$\displaystyle Y_\nu(z)=\frac{-2}{\pi} \int_0 ^\infty \cos(x \cos(t) - \pi \nu/2) \cos(\nu t) \mathrm d t$

and, in particular,

$\displaystyle J_0(x)=\frac{2}{\pi} \int_0^\infty \sin(x \cosh(t)) \mathrm d t$

$\displaystyle Y_0(x)=\frac{-2}{\pi} \int_0^\infty \cos(x \cosh(t)) \mathrm d t$

Also,

$\displaystyle J_\nu(x)=\frac{2 (x/2)^{-\nu}}{\pi^{1/2} \Gamma(1/2-\nu)}

\int_1^\infty \frac{\sin(xt)~ \mathrm d t} { (t^2-1)^{\nu+1/2}}$

$\displaystyle Y_\nu(x)=-\frac{2 (x/2)^{-\nu}}{\pi^{1/2} \Gamma(1/2-\nu)}

\int_1^\infty \frac{\cos(xt)~ \mathrm d t} { (t^2-1)^{\nu+1/2}}$

However the reason of the suggested restriction $x>0$ is not clear. Peerhaps, these expressions can be used to deduce the expansion suitable for the numerical implementation.

Expansion of $J_\nu$ at zero

$\!\!\!\!\!\!\!\!\!\!\!\!\! (9.1.10) ~ ~ ~ \displaystyle

J_\nu(z)=\left(\frac{z}{2}\right)^{\!\nu}~ \sum_{k=0}^{\infty} ~ \frac{(-z^2/4)^k} {k!~ (\nu\!+\!k)!}$

Expansion of $Y_n$ at zero

The similar expansion for $Y_n$ at natural $n\!+\!1$ looks ugly:

$\!\!\!\!\!\!\!\!\!\!\!\!\! (\mathrm{GR} 8.403) ~ ~ ~ \displaystyle

\pi Y_n(z)= 2 J_n(z) ( \ln(z/2) + C ) - \sum_{k=0}^{n-1} \frac{(n\!-\!k\!-\!1)!}{k!} (z/2)^{2k-n} -$

$ \displaystyle

- (z/2)^n \frac{1}{n!} \sum_{k=1}^n \frac{1}{k} - \sum_{k=0}^{\infty} \frac{(-1)^k (z/2)^{n+2k}}{k! ~(n\!+\!k)!} \left( \sum_{m=1}^{n+k} \frac{1}{m} + \sum_{m=1}^k \frac{1}{m} \right) $ where $C$ is Euler constant, called also EulerGamma

$\displaystyle

C=- \int_0^\infty \exp(-t)~ \ln(t) ~\mathrm d t \approx 0.57721566490$

Up to year 2012, no beautiful representation for the expansion coefficients is available.

Expansion at infinity by Gradshtein,Ryzhik

Gradshtein,Ryzhik suggest the following expansions (See 8.4.5.5)

$ \displaystyle J_{\pm \nu}(z)=$

$ \displaystyle = \sqrt{\frac{2}{\pi z}} \cos\!\Big(z-(\pm 2 \nu\!+\!1)\pi/4 \Big) ~ \left( ~ \sum_{k=0}^{n-1} ~

\left(\frac{-1}{4z^2}\right)^{\!\! k} \frac{ \Gamma(\nu+2k+1/2)}{(2k)!~ \Gamma(\nu-2k+1/2)} + R_1 \right) -$

$ \displaystyle - \sqrt{\frac{2}{\pi z}} \sin\!\Big(z-(\pm 2 \nu\!+\!1)\pi/4 \Big)~ \left( \frac{1}{2z} ~ \sum_{k=0}^{n-1} ~

\left(\frac{-1}{4z^2}\right)^{\!\! k} \frac{ \Gamma(\nu+2k+3/2)}{(2k\!+\!1)! ~\Gamma(\nu-2k-1/2)} + R_2 \right)$

$ \displaystyle Y_{\pm \nu}(z)=$

$ \displaystyle = \sqrt{\frac{2}{\pi z}} \sin\!\Big(z-(\pm 2 \nu\!+\!1)\pi/4 \Big) ~ \left( ~ \sum_{k=0}^{n-1} ~

\left(\frac{-1}{4z^2}\right)^{\!\! k} \frac{ \Gamma(\nu+2k+1/2)}{(2k)!~ \Gamma(\nu-2k+1/2)} + R_1 \right) +$

$ \displaystyle + \sqrt{\frac{2}{\pi z}} \cos\!\Big(z-(\pm 2 \nu\!+\!1)\pi/4 \Big)~ \left( \frac{1}{2z} ~ \sum_{k=0}^{n-1} ~

\left(\frac{-1}{4z^2}\right)^{\!\! k} \frac{ \Gamma(\nu+2k+3/2)}{(2k\!+\!1)! ~\Gamma(\nu-2k-1/2)} + R_2 \right)$

$ \displaystyle H_{\nu}(z)=$

$ \displaystyle = \sqrt{\frac{2}{\pi z}} \exp\!\Big(z-(\pm 2 \nu\!+\!1)\pi/4 \Big) ~ \left( ~ \sum_{k=0}^{n-1} ~

\left(\frac{\mathrm i}{2z}\right)^{\!\! k} \frac{ \Gamma(\nu+n+1/2)}{(2n)!~ \Gamma(\nu-n+1/2)} + \theta_1 \left(\frac{\mathrm i}{2z}\right)^{\!\! k} \frac{ \Gamma(\nu+n+1/2)}{(2n)!~ \Gamma(\nu-n+1/2)} \right) $

$ \displaystyle

|R_1|< \left| \frac{\Gamma(\nu+2n+1/2)}{(2z)^{2n} ~ (2n)! ~ \Gamma(\nu-2n+1/2)} \right|$

$ \displaystyle

|R_2|< \left| \frac{\Gamma(\nu+2n+3/2)}{(2z)^{2n+1} ~ (2n\!+\!1)! ~ \Gamma(\nu-2n-1/2)} \right|$

while $\Im(z)\le 0$, the esitmate $|\theta_1| < 1 $

For half–natural $\nu$, the singularity of $\Gamma$ terminates the series and they become finite sums.

Expansion at infinity from Abramowitz,Stegun

Let $\mu=4 \nu^2$. Define two series $P_\nu(z)$ and $Q_\nu(z)$ with

$\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\! (9.2.09) ~ ~ ~ \displaystyle

P_\nu(z)=1 -\frac{(\mu\!-\!1)(\mu\!-\!9)}{2! ~ (8z)^2} +\frac{(\mu\!-\!1)(\mu\!-\!9)(\mu\!-\!25)(\mu\!-\!49)}{4! ~ (8z)^4} - ..$

$\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\! (9.2.10) ~ ~ ~ \displaystyle

Q_\nu(z)= \frac{(\mu\!-\!1)}{1! ~ (8z)} +\frac{(\mu\!-\!1)(\mu\!-\!9)(\mu\!-\!25)}{3! ~ (8z)^3} - \frac{(\mu\!-\!1)(\mu\!-\!9)(\mu\!-\!25)(\mu\!-\!49)(\mu\!-\!81)}{5! ~ (8z)^5} +..$ Let $x=z-\Big(\nu/2+\pi/4\Big)\pi~$; then

$\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\! (9.2.05) ~ ~ ~ \displaystyle

J_\nu(z_)=\sqrt{\frac{2}{\pi z}}\Big(P_\nu(z) \cos(z)-Q_\nu(z) \sin(z) \Big)$

$\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\! (9.2.06) ~ ~ ~ \displaystyle

Y_\nu(z_)=\sqrt{\frac{2}{\pi z}}\Big(P_\nu(z) \sin(z)+Q_\nu(z) \cos(z) \Big)$

$\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\! (9.2.07) ~ ~ ~ \displaystyle

H_\nu(z_)=\sqrt{\frac{2}{\pi z}}\Big(P_\nu(z) + \mathrm i ~ Q_\nu(z) \Big) \mathrm e^{\mathrm i z}$

while $|\mathrm{Arg}(z)|<\pi$.

Asymptotic expansion for modulus and phase

$ \!\!\!\!\!\!\!\!\! (9.2.19) ~ ~ ~ J_\nu(z)= M \cos(\theta) ~~, ~~$ $~ Y_\nu(z)= M \sin(\theta) ~~, ~~$

$ \!\!\!\!\!\!\!\!\! (9.2.17) ~ ~ ~ M = \sqrt{J_\nu(z)^2 + Y_\nu(z)^2}$

In certain range, while $ |\Re(\theta)|<\pi$, also

$ \!\!\!\!\!\!\!\!\! (9.2.0) ~ ~ ~ \theta = \mathrm{atan2}( Y_\nu(z)/ J_\nu(z))$

$M$ is called "modulus" and $\theta$ is called "argument" . Let $\mu=4\nu^2$.

At large values of the argument, it worth to expand $M$ and $\theta$ instead of $J$:

$ \!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\! (9.2.28) \displaystyle ~ ~ ~ M^2 = \frac{2}{\pi z} \left( 1 +

\frac{1}{2} \frac{\mu \!-\!1}{(2z)^2} + \frac{1\cdot 3}{2 \cdot 4} \frac{(\mu \!-\!1)(\mu\!-\!9)}{(2z)^4} + \frac{1\cdot 3 \cdot 5}{2 \cdot 4 \cdot 6} \frac{(\mu\!-\!1)(\mu\!-\!9)(\mu\!-\!25)}{(2z)^6} +.. \right) $

$ \!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\! (9.2.29) \displaystyle ~ ~ ~ \theta = z - \left(\frac{\nu}{2}+\frac{1}{4}\right) \pi

+ \frac{\mu\!-\!1}{2(2z)} + \frac{ (\mu\!-\!1) (\mu-25)}{6(4z)^3} + \frac{ (\mu\!-\!1) (\mu^2-114\mu+1073)}{5(4z)^5} $ $+$ $~\, \displaystyle \frac{ (\mu\!-\!1) (5\mu^3-1535\mu^2+54703\mu-375733) }{14(4z)^7} + ..$

Integrals

Some examples with Mathematica are shown below:

Integrate[BesselJ[0, x y] y, {y, 0, 1}] $~$ does formula

$\displaystyle \int_0^1 J_0(x y) \, y\, \mathrm d y = \frac{ J_1(x)}{x}$

Simplify[Integrate[1/Sqrt[1-x^2] BesselJ[0, k x], {x, 0, 1}], k > 0] $~$ does formula

$\displaystyle \int_0^1 \frac{1}{\sqrt{1-x^2}}\, J_0(kx) \, \mathrm d x = \frac{\pi}{2} J_0(k/2)^2$

Simplify[Integrate[1/Sqrt[1-x^2] BesselJ[0, k x] x, {x, 0, 1}], k > 0] $~$ does formula

$\displaystyle \int_0^1 \frac{1}{\sqrt{1-x^2}} \, J_0(kx) \, x \, \mathrm d x = \frac{\sin(k)}{k}$

These expressions can be used for testing of the numerical implementations of the Bessel transform. Wolfram suggests more functions such that their Bessel transforms are expressed in terms of special functions at http://mathworld.wolfram.com/HankelTransform.html

Derivatives of the Bessel

D[BesselJ[v, x], x] $~ ~$ gives

1/2 (BesselJ[-1 + v, x] - BesselJ[1 + v, x])

In "short" notations, this I can be written with

${J_v}^{\prime}(x)=\frac{1}{2}\big( J_{v-1}(x)-J_{v+1}(x) \big)$

Zeros of the Bessel

Zeros of the Bessel function are demote with identifier BesselJZero; first argument indicates the order of the Bessel function; second one base sense of the radial coordinate. For the numerical evaluation, the asymptotic expression can be used,

zo[k_] = Series[BesselJZero[v, k], {k, Infinity, 1}]

$\pi \left(k+\frac{v}{2}-\frac{1}{4}\right)-\frac{(2

  v-1) (2 v+1)}{8 \pi 
  \left(k+\frac{v}{2}-\frac{1}{4}\right)}$

For the precise evaluation at integer $v$, three steps of the Newton iteration to adjust the value seem to be sufficient, at least for the "double" precision. The C++ implementation can be denoted jnz or, as in Mathematica, BesselJZero. The example and the test:

#include<math.h>
#include<stdio.h>
#define DB double
#define DO(x,y) for(x=0;x<y;x++)

DB jnp(int n,DB x){ return .5*( jn(n-1,x)-jn(n+1,x) ) ; } // Derivative of n th Bessel

DB jnz(int v, int k){ DB x,t; t=M_PI*(k+.5*v-.25); x= t - (v*v-.25)*.5/t;
                x-= jn(v,x)/jnp(v,x); // Newton adjustment of the root
                x-= jn(v,x)/jnp(v,x);
                x-= jn(v,x)/jnp(v,x);
                return x; } // the k th zero of v th Bessel

int main(){ int n,v,k; DB x,y; // evaluation of zeros of Bessel of order 0,1,2,3 and the testing
for(n=1;n<101;n++){ printf("%3d",n);
                DO(v,4){x=jnz(v,n); y=jn(v,x);
                        printf("%21.16lf %19.16lf",x,y);
                        }
                printf("\n");
                }
return 0;
}

The first 10 lines of the output are copupasted below:

 1   2.4048255576957729 -0.0000000000000001   3.8317059702075125 -0.0000000000000000   5.1356223018406828 -0.0000000000000001   6.3801618959239841 -0.0000000000000002
 2   5.5200781102863106 -0.0000000000000000   7.0155866698156188  0.0000000000000000   8.4172441403998643 -0.0000000000000001   9.7610231299816697 -0.0000000000000000
 3   8.6537279129110125 -0.0000000000000001  10.1734681350627216  0.0000000000000001  11.6198411721490587  0.0000000000000002  13.0152007216984344 -0.0000000000000000
 4  11.7915344390142813 -0.0000000000000001  13.3236919363142228 -0.0000000000000001  14.7959517823512599 -0.0000000000000002  16.2234661603187682  0.0000000000000000
 5  14.9309177084877867 -0.0000000000000001  16.4706300508776344 -0.0000000000000003  17.9598194949878263  0.0000000000000000  19.4094152264350122 -0.0000000000000001
 6  18.0710639679109235  0.0000000000000002  19.6158585104682430  0.0000000000000002  21.1169970530218443 -0.0000000000000002  22.5827295931044425  0.0000000000000001
 7  21.2116366298792585  0.0000000000000001  22.7600843805927724 -0.0000000000000001  24.2701123135731009  0.0000000000000003  25.7481666992949769  0.0000000000000001
 8  24.3524715307493018 -0.0000000000000001  25.9036720876183821 -0.0000000000000001  27.4205735499845566 -0.0000000000000001  28.9083507809217579  0.0000000000000000
 9  27.4934791320402532  0.0000000000000002  29.0468285349168553 -0.0000000000000000  30.5692044955163986 -0.0000000000000002  32.0648524070977103 -0.0000000000000001
10  30.6346064684319757  0.0000000000000001  32.1896799109744052  0.0000000000000002  33.7165195092226995 -0.0000000000000001  35.2186707386101148  0.0000000000000000

References

http://en.wikipedia.org/wiki/Bessel_function

http://en.citizendium.org/wiki/Bessel_functions

Keywords

BesselJ0, BesselJ1, BesselY0, BesselY1,