WSJT-X/boost/libs/math/doc/sf/bessel_ik.qbk

[section:mbessel Modified Bessel Functions of the First and Second Kinds]

[h4 Synopsis]

`#include <boost/math/special_functions/bessel.hpp>`

   template <class T1, class T2>
   ``__sf_result`` cyl_bessel_i(T1 v, T2 x);

   template <class T1, class T2, class ``__Policy``>
   ``__sf_result`` cyl_bessel_i(T1 v, T2 x, const ``__Policy``&);

   template <class T1, class T2>
   ``__sf_result`` cyl_bessel_k(T1 v, T2 x);
   
   template <class T1, class T2, class ``__Policy``>
   ``__sf_result`` cyl_bessel_k(T1 v, T2 x, const ``__Policy``&);
   
   
[h4 Description]

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_k(v, x) = K[sub v](x)

where:

[equation mbessel2]

[equation mbessel3]

The return type of these functions is computed using the __arg_promotion_rules
when T1 and T2 are different types.  The functions are also optimised for the
relatively common case that T1 is an integer.

[optional_policy]

The functions return the result of __domain_error whenever the result is
undefined or complex.  For __cyl_bessel_j this occurs when `x < 0` and v is not
an integer, or when `x == 0` and `v != 0`.  For __cyl_neumann this occurs
when `x <= 0`.

The following graph illustrates the exponential behaviour of I[sub v].

[graph cyl_bessel_i]

The following graph illustrates the exponential decay of K[sub v].

[graph cyl_bessel_k]

[h4 Testing]

There are two sets of test values: spot values calculated using
[@http://functions.wolfram.com functions.wolfram.com],
and a much larger set of tests computed using
a simplified version of this implementation
(with all the special case handling removed).

[h4 Accuracy]

The following tables show how the accuracy of these functions
varies on various platforms, along with comparison to other libraries.  
Note that only results for the widest floating-point type on the 
system are given, as narrower types have __zero_error.  All values
are relative errors in units of epsilon.  Note that our test suite
includes some fairly extreme inputs which results in most of the worst
problem cases in other libraries:

[table_cyl_bessel_i_integer_orders_]

[table_cyl_bessel_i]

[table_cyl_bessel_k_integer_orders_]

[table_cyl_bessel_k]

[h4 Implementation]

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
['x > 0].

For I[sub v][space] with v equal to 0, 1 or 0.5 are handled as special cases.

The 0 and 1 cases use minimax rational approximations on
finite and infinite intervals. The coefficients are from:

* J.M. Blair and C.A. Edwards, ['Stable rational minimax approximations
    to the modified Bessel functions I_0(x) and I_1(x)], Atomic Energy of Canada
    Limited Report 4928, Chalk River, 1974.
* S. Moshier, ['Methods and Programs for Mathematical Functions],
    Ellis Horwood Ltd, Chichester, 1989.    

While the 0.5 case is a simple trigonometric function:

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:

[equation mbessel5]

starting from K[sub 0][space] and K[sub 1][space] 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:

[equation mbessel17]

In the general case, we first normalize [nu][space] 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
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), 
which needs two continued fractions as well as the Wronskian:

[equation mbessel11]

[equation mbessel12]

[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)). 
Their convergence rates depend on ['x], therefore we need
different strategies for large ['x] and small ['x].

['x > v], CF1 needs O(['x]) iterations to converge, CF2 converges rapidly.

['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]
can be calculated by

[equation mbessel13]

where

[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
    of real order and complex argument to selected accuracy], Computer Physics
    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 mbessel15]

where

[equation mbessel16]

f[sub k][space] and h[sub k][space]
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).
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]

[/ 
  Copyright 2006 John Maddock, Paul A. Bristow and Xiaogang Zhang.
  Distributed under the Boost Software License, Version 1.0.
  (See accompanying file LICENSE_1_0.txt or copy at
  http://www.boost.org/LICENSE_1_0.txt).
]