[section:mbessel Modified Bessel Functions of the First and Second Kinds] [h4 Synopsis] template ``__sf_result`` cyl_bessel_i(T1 v, T2 x); template ``__sf_result`` cyl_bessel_i(T1 v, T2 x, const ``__Policy``&); template ``__sf_result`` cyl_bessel_k(T1 v, T2 x); template ``__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_pomotion_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 a comparison to the __gsl library. 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. [table Errors Rates in cyl_bessel_i [[Significand Size] [Platform and Compiler] [I[sub v]] ] [[53] [Win32 / Visual C++ 8.0] [Peak=10 Mean=3.4 GSL Peak=6000] ] [[64] [Red Hat Linux IA64 / G++ 3.4] [Peak=11 Mean=3] ] [[64] [SUSE Linux AMD64 / G++ 4.1] [Peak=11 Mean=4] ] [[113] [HP-UX / HP aCC 6] [Peak=15 Mean=4] ] ] [table Errors Rates in cyl_bessel_k [[Significand Size] [Platform and Compiler] [K[sub v]] ] [[53] [Win32 / Visual C++ 8.0] [Peak=9 Mean=2 GSL Peak=9] ] [[64] [Red Hat Linux IA64 / G++ 3.4] [Peak=10 Mean=2] ] [[64] [SUSE Linux AMD64 / G++ 4.1] [Peak=10 Mean=2] ] [[113] [HP-UX / HP aCC 6] [Peak=12 Mean=5] ] ] [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). ]