Difference between revisions of "Bessel function"
Jump to navigation
Jump to search
m (Text replacement - "\$([^\$]+)\$" to "\\(\1\\)") |
|||
| Line 1: | Line 1: | ||
| − | '''Bessel function''' referes to a solution |
+ | '''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: |
with the following asymptotic behaviour: |
||
| − | + | \(\displaystyle |
|
| − | f(x) \approx x^\nu \left( \frac{2^{-\nu}}{\mathrm{Factorial}(\nu)}+ O(x^2) \right) |
+ | f(x) \approx x^\nu \left( \frac{2^{-\nu}}{\mathrm{Factorial}(\nu)}+ O(x^2) \right)\) |
In particular, |
In particular, |
||
| − | For |
+ | For \(\nu=0\), the initial condition is assumed \(~f(0)=1~\), \(~f'(0)=0\) , and |
| − | For |
+ | For \(\nu=1\), the initial condition is assumed \(~f(0)=0~\), \(~f'(0)=1\); |
| − | Such a solution |
+ | Such a solution \(f\) denoted with \(\mathrm{BesselJ}[\nu,x]\) or with \(J_\nu(x)\). |
==BesselJ== |
==BesselJ== |
||
| − | Due to singularity of the equation at |
+ | 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]] |
Many formulas about the Bessel functions below are borrowed from the handbook by [[Abramowirtz,Stegun]] |
||
<ref name="a"> |
<ref name="a"> |
||
| Line 28: | Line 28: | ||
==Integral representations== |
==Integral representations== |
||
| − | : |
+ | : \( \!\!\!\!\!\!\!\!\!\! (9.1.20) ~ ~ ~ \displaystyle |
J_\nu(z) = \frac{(z/2)^{\nu}}{\pi^{1/2} ~(\nu-1/2)!} |
J_\nu(z) = \frac{(z/2)^{\nu}}{\pi^{1/2} ~(\nu-1/2)!} |
||
~ |
~ |
||
| Line 34: | Line 34: | ||
~ |
~ |
||
\cos(z \cos(t)) \sin(t)^{2 \nu} ~t~ \mathrm d t |
\cos(z \cos(t)) \sin(t)^{2 \nu} ~t~ \mathrm d t |
||
| + | \) |
||
| − | $ |
||
===Sonin representation=== |
===Sonin representation=== |
||
| Line 42: | Line 42: | ||
</ref> suggest that |
</ref> 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, |
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, |
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}} |
+ | \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}} |
+ | \int_1^\infty \frac{\cos(xt)~ \mathrm d t} { (t^2-1)^{\nu+1/2}}\) |
| − | However the reason of the suggested restriction |
+ | 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. |
Peerhaps, these expressions can be used to deduce the expansion suitable for the numerical implementation. |
||
| − | ==Expansion of |
+ | ==Expansion of \(J_\nu\) at zero== |
| − | : |
+ | : \(\!\!\!\!\!\!\!\!\!\!\!\!\! (9.1.10) ~ ~ ~ \displaystyle |
J_\nu(z)=\left(\frac{z}{2}\right)^{\!\nu}~ |
J_\nu(z)=\left(\frac{z}{2}\right)^{\!\nu}~ |
||
\sum_{k=0}^{\infty} ~ |
\sum_{k=0}^{\infty} ~ |
||
\frac{(-z^2/4)^k} |
\frac{(-z^2/4)^k} |
||
| − | {k!~ (\nu\!+\!k)!} |
+ | {k!~ (\nu\!+\!k)!}\) |
| − | ==Expansion of |
+ | ==Expansion of \(Y_n\) at zero== |
| − | The similar expansion for |
+ | 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} - |
+ | \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)!} |
- (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( |
\left( |
||
\sum_{m=1}^{n+k} \frac{1}{m} + \sum_{m=1}^k \frac{1}{m} \right) |
\sum_{m=1}^{n+k} \frac{1}{m} + \sum_{m=1}^k \frac{1}{m} \right) |
||
| + | \) |
||
| − | $ |
||
| − | where |
+ | where \(C\) is [[Euler constant]], called also [[EulerGamma]] |
| − | : |
+ | : \(\displaystyle |
| − | C=- \int_0^\infty \exp(-t)~ \ln(t) ~\mathrm d t \approx 0.57721566490 |
+ | 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. |
Up to year 2012, no beautiful representation for the expansion coefficients is available. |
||
| Line 98: | Line 98: | ||
!--> |
!--> |
||
| − | : |
+ | : \( \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} |
\left(\frac{-1}{4z^2}\right)^{\!\! k} |
||
| − | \frac{ \Gamma(\nu+2k+1/2)}{(2k)!~ \Gamma(\nu-2k+1/2)} + R_1 \right) - |
+ | \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} |
\left(\frac{-1}{4z^2}\right)^{\!\! k} |
||
| − | \frac{ \Gamma(\nu+2k+3/2)}{(2k\!+\!1)! ~\Gamma(\nu-2k-1/2)} + R_2 \right) |
+ | \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} |
\left(\frac{-1}{4z^2}\right)^{\!\! k} |
||
| − | \frac{ \Gamma(\nu+2k+1/2)}{(2k)!~ \Gamma(\nu-2k+1/2)} + R_1 \right) + |
+ | \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} |
\left(\frac{-1}{4z^2}\right)^{\!\! k} |
||
| − | \frac{ \Gamma(\nu+2k+3/2)}{(2k\!+\!1)! ~\Gamma(\nu-2k-1/2)} + R_2 \right) |
+ | \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} |
\left(\frac{\mathrm i}{2z}\right)^{\!\! k} |
||
\frac{ \Gamma(\nu+n+1/2)}{(2n)!~ \Gamma(\nu-n+1/2)} + |
\frac{ \Gamma(\nu+n+1/2)}{(2n)!~ \Gamma(\nu-n+1/2)} + |
||
\theta_1 \left(\frac{\mathrm i}{2z}\right)^{\!\! k} |
\theta_1 \left(\frac{\mathrm i}{2z}\right)^{\!\! k} |
||
| − | \frac{ \Gamma(\nu+n+1/2)}{(2n)!~ \Gamma(\nu-n+1/2)} \right) |
+ | \frac{ \Gamma(\nu+n+1/2)}{(2n)!~ \Gamma(\nu-n+1/2)} \right) \) |
| − | : |
+ | : \( \displaystyle |
|R_1|< \left| |
|R_1|< \left| |
||
\frac{\Gamma(\nu+2n+1/2)}{(2z)^{2n} ~ (2n)! ~ \Gamma(\nu-2n+1/2)} |
\frac{\Gamma(\nu+2n+1/2)}{(2z)^{2n} ~ (2n)! ~ \Gamma(\nu-2n+1/2)} |
||
| − | \right| |
+ | \right|\) |
| − | : |
+ | : \( \displaystyle |
|R_2|< \left| |
|R_2|< \left| |
||
\frac{\Gamma(\nu+2n+3/2)}{(2z)^{2n+1} ~ (2n\!+\!1)! ~ \Gamma(\nu-2n-1/2)} |
\frac{\Gamma(\nu+2n+3/2)}{(2z)^{2n+1} ~ (2n\!+\!1)! ~ \Gamma(\nu-2n-1/2)} |
||
| − | \right| |
+ | \right|\) |
| − | : while |
+ | : while \(\Im(z)\le 0\), the esitmate \(|\theta_1| < 1 \) |
| − | For half–natural |
+ | For half–natural \(\nu\), the singularity of \(\Gamma\) terminates the series and they become finite sums. |
==Expansion at infinity from [[Abramowitz,Stegun]]== |
==Expansion at infinity from [[Abramowitz,Stegun]]== |
||
| − | Let |
+ | Let \(\mu=4 \nu^2\). Define two series \(P_\nu(z)\) and \(Q_\nu(z)\) with |
| − | : |
+ | : \(\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\! (9.2.09) ~ ~ ~ \displaystyle |
P_\nu(z)=1 |
P_\nu(z)=1 |
||
-\frac{(\mu\!-\!1)(\mu\!-\!9)}{2! ~ (8z)^2} |
-\frac{(\mu\!-\!1)(\mu\!-\!9)}{2! ~ (8z)^2} |
||
+\frac{(\mu\!-\!1)(\mu\!-\!9)(\mu\!-\!25)(\mu\!-\!49)}{4! ~ (8z)^4} |
+\frac{(\mu\!-\!1)(\mu\!-\!9)(\mu\!-\!25)(\mu\!-\!49)}{4! ~ (8z)^4} |
||
| − | - .. |
+ | - ..\) |
| − | : |
+ | : \(\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\! (9.2.10) ~ ~ ~ \displaystyle |
Q_\nu(z)= |
Q_\nu(z)= |
||
\frac{(\mu\!-\!1)}{1! ~ (8z)} |
\frac{(\mu\!-\!1)}{1! ~ (8z)} |
||
+\frac{(\mu\!-\!1)(\mu\!-\!9)(\mu\!-\!25)}{3! ~ (8z)^3} |
+\frac{(\mu\!-\!1)(\mu\!-\!9)(\mu\!-\!25)}{3! ~ (8z)^3} |
||
- \frac{(\mu\!-\!1)(\mu\!-\!9)(\mu\!-\!25)(\mu\!-\!49)(\mu\!-\!81)}{5! ~ (8z)^5} |
- \frac{(\mu\!-\!1)(\mu\!-\!9)(\mu\!-\!25)(\mu\!-\!49)(\mu\!-\!81)}{5! ~ (8z)^5} |
||
| − | +.. |
+ | +..\) |
| − | Let |
+ | Let \(x=z-\Big(\nu/2+\pi/4\Big)\pi~\); |
then |
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) |
+ | 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) |
+ | 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} |
+ | H_\nu(z_)=\sqrt{\frac{2}{\pi z}}\Big(P_\nu(z) + \mathrm i ~ Q_\nu(z) \Big) \mathrm e^{\mathrm i z}\) |
| − | while |
+ | while \(|\mathrm{Arg}(z)|<\pi\). |
==Asymptotic expansion for modulus and phase== |
==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 |
+ | 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" |
|
<ref name="a"> |
<ref name="a"> |
||
http://people.math.sfu.ca/~cbm/aands/page_365.htm |
http://people.math.sfu.ca/~cbm/aands/page_365.htm |
||
Abramovitz, Stegun. Handbook on mathematical functions. |
Abramovitz, Stegun. Handbook on mathematical functions. |
||
| − | </ref>. Let |
+ | </ref>. Let \(\mu=4\nu^2\). |
| − | At large values of the argument, it worth to expand |
+ | 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}{2} \frac{\mu \!-\!1}{(2z)^2} + |
||
\frac{1\cdot 3}{2 \cdot 4} \frac{(\mu \!-\!1)(\mu\!-\!9)}{(2z)^4} + |
\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} +.. |
\frac{1\cdot 3 \cdot 5}{2 \cdot 4 \cdot 6} \frac{(\mu\!-\!1)(\mu\!-\!9)(\mu\!-\!25)}{(2z)^6} +.. |
||
\right) |
\right) |
||
| + | \) |
||
| − | $ |
||
| − | : |
+ | : \( \!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\! (9.2.29) \displaystyle ~ ~ ~ \theta = z - \left(\frac{\nu}{2}+\frac{1}{4}\right) \pi |
+ \frac{\mu\!-\!1}{2(2z)} |
+ \frac{\mu\!-\!1}{2(2z)} |
||
+ \frac{ (\mu\!-\!1) (\mu-25)}{6(4z)^3} |
+ \frac{ (\mu\!-\!1) (\mu-25)}{6(4z)^3} |
||
| − | + \frac{ (\mu\!-\!1) (\mu^2-114\mu+1073)}{5(4z)^5} |
+ | + \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== |
==Integrals== |
||
| Line 196: | Line 196: | ||
<nomathjax> |
<nomathjax> |
||
Integrate[BesselJ[0, x y] y, {y, 0, 1}] |
Integrate[BesselJ[0, x y] y, {y, 0, 1}] |
||
| − | </nomathjax> |
+ | </nomathjax> \(~\) |
does formula |
does formula |
||
| − | + | \(\displaystyle \int_0^1 J_0(x y) \, y\, \mathrm d y = \frac{ J_1(x)}{x}\) |
|
<nomathjax> |
<nomathjax> |
||
Simplify[Integrate[1/Sqrt[1-x^2] BesselJ[0, k x], {x, 0, 1}], k > 0] |
Simplify[Integrate[1/Sqrt[1-x^2] BesselJ[0, k x], {x, 0, 1}], k > 0] |
||
| − | </nomathjax> |
+ | </nomathjax> \(~\) |
does formula |
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\) |
|
<nomathjax> |
<nomathjax> |
||
Simplify[Integrate[1/Sqrt[1-x^2] BesselJ[0, k x] x, {x, 0, 1}], k > 0] |
Simplify[Integrate[1/Sqrt[1-x^2] BesselJ[0, k x] x, {x, 0, 1}], k > 0] |
||
| − | </nomathjax> |
+ | </nomathjax> \(~\) |
does formula |
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]]. |
These expressions can be used for testing of the numerical implementations of the [[Bessel transform]]. |
||
| Line 220: | Line 220: | ||
==Derivatives of the Bessel== |
==Derivatives of the Bessel== |
||
| − | D[BesselJ[v, x], x] |
+ | D[BesselJ[v, x], x] \(~ ~\) gives |
1/2 (BesselJ[-1 + v, x] - BesselJ[1 + v, x]) |
1/2 (BesselJ[-1 + v, x] - BesselJ[1 + v, x]) |
||
| Line 226: | Line 226: | ||
In "short" notations, this I can be written with |
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) |
J_{v-1}(x)-J_{v+1}(x) |
||
| − | \big) |
+ | \big)\) |
==Zeros of the Bessel== |
==Zeros of the Bessel== |
||
| Line 237: | Line 237: | ||
zo[k_] = Series[BesselJZero[v, k], {k, Infinity, 1}] |
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 |
v-1) (2 v+1)}{8 \pi |
||
| − | \left(k+\frac{v}{2}-\frac{1}{4}\right)} |
+ | \left(k+\frac{v}{2}-\frac{1}{4}\right)}\) |
| − | For the precise evaluation at integer |
+ | 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: |
The [[C++]] implementation can be denoted [[jnz]] or, as in Mathematica, [[BesselJZero]]. The example and the test: |
||
<poem><nomathjax><nowiki> |
<poem><nomathjax><nowiki> |
||
Latest revision as of 18:25, 30 July 2019
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