-
[section:mbessel Modified Bessel Functions of the First and Second Kinds]
[h4 Synopsis]
The functions __cyl_bessel_i and __cyl_bessel_k return the result of the
modified Bessel functions of the first and second kind respectively:
-cyl_bessel_i(v, x) = I[sub v](x)
+[:cyl_bessel_i(v, x) = I[sub v](x)]
-cyl_bessel_k(v, x) = K[sub v](x)
+[:cyl_bessel_k(v, x) = K[sub v](x)]
where:
The following error plot are based on an exhaustive search of the functions domain for I0, I1, K0, and K1,
MSVC-15.5 at `double` precision, and GCC-7.1/Ubuntu for `long double` and `__float128`.
-
[graph i0__double]
[graph i0__80_bit_long_double]
The following are handled as special cases first:
-When computing I[sub v][space] for ['x < 0], then [nu][space] must be an integer
-or a domain error occurs. If [nu][space] is an integer, then the function is
-odd if [nu][space] is odd and even if [nu][space] is even, and we can reflect to
+When computing I[sub v] for ['x < 0], then [nu] must be an integer
+or a domain error occurs. If [nu] is an integer, then the function is
+odd if [nu] is odd and even if [nu] is even, and we can reflect to
['x > 0].
-For I[sub v][space] with v equal to 0, 1 or 0.5 are handled as special cases.
+For I[sub v] with v equal to 0, 1 or 0.5 are handled as special cases.
The 0 and 1 cases use polynomial approximations on
finite and infinite intervals. The approximating forms
The 0.5 case is a simple trigonometric function:
-I[sub 0.5](x) = sqrt(2 / [pi]x) * sinh(x)
+[:I[sub 0.5](x) = sqrt(2 / [pi]x) * sinh(x)]
-For K[sub v][space] with /v/ an integer, the result is calculated using the
-recurrence relation:
+For K[sub v] with /v/ an integer, the result is calculated using the recurrence relation:
[equation mbessel5]
-starting from K[sub 0][space] and K[sub 1][space] which are calculated
+starting from K[sub 0] and K[sub 1] which are calculated
using rational the approximations above. These rational approximations are
accurate to around 19 digits, and are therefore only used when T has
no more than 64 binary digits of precision.
-When /x/ is small compared to /v/, I[sub v]x[space] is best computed directly from the series:
+When /x/ is small compared to /v/, I[sub v]x is best computed directly from the series:
[equation mbessel17]
-In the general case, we first normalize [nu][space] to \[[^0, [inf]])
+In the general case, we first normalize [nu] to \[[^0, [inf]])
with the help of the reflection formulae:
[equation mbessel9]
[equation mbessel10]
-Let [mu][space] = [nu] - floor([nu] + 1/2), then [mu][space] is the fractional part of
-[nu][space] such that |[mu]| <= 1/2 (we need this for convergence later). The idea is to
+Let [mu] = [nu] - floor([nu] + 1/2), then [mu] is the fractional part of
+[nu] such that |[mu]| <= 1/2 (we need this for convergence later). The idea is to
calculate K[sub [mu]](x) and K[sub [mu]+1](x), and use them to obtain
I[sub [nu]](x) and K[sub [nu]](x).
The algorithm is proposed by Temme in
-N.M. Temme, ['On the numerical evaluation of the modified bessel function
- of the third kind], Journal of Computational Physics, vol 19, 324 (1975),
+[:N.M. Temme, ['On the numerical evaluation of the modified bessel function
+ of the third kind], Journal of Computational Physics, vol 19, 324 (1975),]
which needs two continued fractions as well as the Wronskian:
[equation mbessel11]
[equation mbessel8]
The continued fractions are computed using the modified Lentz's method
-(W.J. Lentz, ['Generating Bessel functions in Mie scattering calculations
- using continued fractions], Applied Optics, vol 15, 668 (1976)).
+[:(W.J. Lentz, ['Generating Bessel functions in Mie scattering calculations
+ using continued fractions], Applied Optics, vol 15, 668 (1976)).]
Their convergence rates depend on ['x], therefore we need
different strategies for large ['x] and small ['x].
['x <= v], CF1 converges rapidly, CF2 fails to converge when ['x] [^->] 0.
When ['x] is large (['x] > 2), both continued fractions converge (CF1
-may be slow for really large ['x]). K[sub [mu]][space] and K[sub [mu]+1][space]
+may be slow for really large ['x]). K[sub [mu]] and K[sub [mu]+1]
can be calculated by
[equation mbessel13]
[equation mbessel14]
['S] is also a series that is summed along with CF2, see
-I.J. Thompson and A.R. Barnett, ['Modified Bessel functions I_v and K_v
+[:I.J. Thompson and A.R. Barnett, ['Modified Bessel functions I_v and K_v
of real order and complex argument to selected accuracy], Computer Physics
- Communications, vol 47, 245 (1987).
+ Communications, vol 47, 245 (1987).]
When ['x] is small (['x] <= 2), CF2 convergence may fail (but CF1
works very well). The solution here is Temme's series:
[equation mbessel16]
-f[sub k][space] and h[sub k][space]
+f[sub k] and h[sub k]
are also computed by recursions (involving gamma functions), but the
formulas are a little complicated, readers are referred to
-N.M. Temme, ['On the numerical evaluation of the modified Bessel function
- of the third kind], Journal of Computational Physics, vol 19, 324 (1975).
+[:N.M. Temme, ['On the numerical evaluation of the modified Bessel function
+ of the third kind], Journal of Computational Physics, vol 19, 324 (1975).]
Note: Temme's series converge only for |[mu]| <= 1/2.
K[sub [nu]](x) is then calculated from the forward
recurrence, as is K[sub [nu]+1](x). With these two values and
f[sub [nu]], the Wronskian yields I[sub [nu]](x) directly.
-[endsect]
+[endsect] [/section:mbessel Modified Bessel Functions of the First and Second Kinds]
[/
Copyright 2006 John Maddock, Paul A. Bristow and Xiaogang Zhang.
projects / platform / upstream / boost.git / blobdiff