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 [1]; 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 [2] 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 [3] 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" [1]. 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
- ↑ 1.0 1.1 http://people.math.sfu.ca/~cbm/aands/page_365.htm Abramovitz, Stegun. Handbook on mathematical functions.
- ↑ http://dlmf.nist.gov/10.9 Digital library of mathematical functions
- ↑ http://books.google.co.jp/books?id=aBgFYxKHUjsC&pg=PA859&hl=ja&source=gbs_toc_r&cad=4#v=onepage&q&f=false Izrail Solomonovich Gradshtein, Iosif Moiseevich Ryzhik, Alan Jeffrey, Daniel Zwillinger. Table of Integrals, Series, And Products.
http://en.wikipedia.org/wiki/Bessel_function
http://en.citizendium.org/wiki/Bessel_functions