1[section:mbessel Modified Bessel Functions of the First and Second Kinds] 2 3[h4 Synopsis] 4 5`#include <boost/math/special_functions/bessel.hpp>` 6 7 template <class T1, class T2> 8 ``__sf_result`` cyl_bessel_i(T1 v, T2 x); 9 10 template <class T1, class T2, class ``__Policy``> 11 ``__sf_result`` cyl_bessel_i(T1 v, T2 x, const ``__Policy``&); 12 13 template <class T1, class T2> 14 ``__sf_result`` cyl_bessel_k(T1 v, T2 x); 15 16 template <class T1, class T2, class ``__Policy``> 17 ``__sf_result`` cyl_bessel_k(T1 v, T2 x, const ``__Policy``&); 18 19 20[h4 Description] 21 22The functions __cyl_bessel_i and __cyl_bessel_k return the result of the 23modified Bessel functions of the first and second kind respectively: 24 25[:cyl_bessel_i(v, x) = I[sub v](x)] 26 27[:cyl_bessel_k(v, x) = K[sub v](x)] 28 29where: 30 31[equation mbessel2] 32 33[equation mbessel3] 34 35The return type of these functions is computed using the __arg_promotion_rules 36when T1 and T2 are different types. The functions are also optimised for the 37relatively common case that T1 is an integer. 38 39[optional_policy] 40 41The functions return the result of __domain_error whenever the result is 42undefined or complex. For __cyl_bessel_j this occurs when `x < 0` and v is not 43an integer, or when `x == 0` and `v != 0`. For __cyl_neumann this occurs 44when `x <= 0`. 45 46The following graph illustrates the exponential behaviour of I[sub v]. 47 48[graph cyl_bessel_i] 49 50The following graph illustrates the exponential decay of K[sub v]. 51 52[graph cyl_bessel_k] 53 54[h4 Testing] 55 56There are two sets of test values: spot values calculated using 57[@http://functions.wolfram.com functions.wolfram.com], 58and a much larger set of tests computed using 59a simplified version of this implementation 60(with all the special case handling removed). 61 62[h4 Accuracy] 63 64The following tables show how the accuracy of these functions 65varies on various platforms, along with comparison to other libraries. 66Note that only results for the widest floating-point type on the 67system are given, as narrower types have __zero_error. All values 68are relative errors in units of epsilon. Note that our test suite 69includes some fairly extreme inputs which results in most of the worst 70problem cases in other libraries: 71 72[table_cyl_bessel_i_integer_orders_] 73 74[table_cyl_bessel_i] 75 76[table_cyl_bessel_k_integer_orders_] 77 78[table_cyl_bessel_k] 79 80The following error plot are based on an exhaustive search of the functions domain for I0, I1, K0, and K1, 81MSVC-15.5 at `double` precision, and GCC-7.1/Ubuntu for `long double` and `__float128`. 82 83[graph i0__double] 84 85[graph i0__80_bit_long_double] 86 87[graph i0____float128] 88 89[graph i1__double] 90 91[graph i1__80_bit_long_double] 92 93[graph i1____float128] 94 95[graph k0__double] 96 97[graph k0__80_bit_long_double] 98 99[graph k0____float128] 100 101[graph k1__double] 102 103[graph k1__80_bit_long_double] 104 105[graph k1____float128] 106 107 108[h4 Implementation] 109 110The following are handled as special cases first: 111 112When computing I[sub v] for ['x < 0], then [nu] must be an integer 113or a domain error occurs. If [nu] is an integer, then the function is 114odd if [nu] is odd and even if [nu] is even, and we can reflect to 115['x > 0]. 116 117For I[sub v] with v equal to 0, 1 or 0.5 are handled as special cases. 118 119The 0 and 1 cases use polynomial approximations on 120finite and infinite intervals. The approximating forms 121are based on 122[@http://www.advanpix.com/2015/11/11/rational-approximations-for-the-modified-bessel-function-of-the-first-kind-i0-computations-double-precision/ 123"Rational Approximations for the Modified Bessel Function of the First Kind - I[sub 0](x) for Computations with Double Precision"] 124by Pavel Holoborodko, extended by us to deal with up to 128-bit precision (with different approximations for each target precision). 125 126[equation bessel21] 127 128[equation bessel20] 129 130[equation bessel17] 131 132[equation bessel18] 133 134Similarly we have: 135 136[equation bessel22] 137 138[equation bessel23] 139 140[equation bessel24] 141 142[equation bessel25] 143 144The 0.5 case is a simple trigonometric function: 145 146[:I[sub 0.5](x) = sqrt(2 / [pi]x) * sinh(x)] 147 148For K[sub v] with /v/ an integer, the result is calculated using the recurrence relation: 149 150[equation mbessel5] 151 152starting from K[sub 0] and K[sub 1] which are calculated 153using rational the approximations above. These rational approximations are 154accurate to around 19 digits, and are therefore only used when T has 155no more than 64 binary digits of precision. 156 157When /x/ is small compared to /v/, I[sub v]x is best computed directly from the series: 158 159[equation mbessel17] 160 161In the general case, we first normalize [nu] to \[[^0, [inf]]) 162with the help of the reflection formulae: 163 164[equation mbessel9] 165 166[equation mbessel10] 167 168Let [mu] = [nu] - floor([nu] + 1/2), then [mu] is the fractional part of 169[nu] such that |[mu]| <= 1/2 (we need this for convergence later). The idea is to 170calculate K[sub [mu]](x) and K[sub [mu]+1](x), and use them to obtain 171I[sub [nu]](x) and K[sub [nu]](x). 172 173The algorithm is proposed by Temme in 174[:N.M. Temme, ['On the numerical evaluation of the modified bessel function 175 of the third kind], Journal of Computational Physics, vol 19, 324 (1975),] 176which needs two continued fractions as well as the Wronskian: 177 178[equation mbessel11] 179 180[equation mbessel12] 181 182[equation mbessel8] 183 184The continued fractions are computed using the modified Lentz's method 185[:(W.J. Lentz, ['Generating Bessel functions in Mie scattering calculations 186 using continued fractions], Applied Optics, vol 15, 668 (1976)).] 187Their convergence rates depend on ['x], therefore we need 188different strategies for large ['x] and small ['x]. 189 190['x > v], CF1 needs O(['x]) iterations to converge, CF2 converges rapidly. 191 192['x <= v], CF1 converges rapidly, CF2 fails to converge when ['x] [^->] 0. 193 194When ['x] is large (['x] > 2), both continued fractions converge (CF1 195may be slow for really large ['x]). K[sub [mu]] and K[sub [mu]+1] 196can be calculated by 197 198[equation mbessel13] 199 200where 201 202[equation mbessel14] 203 204['S] is also a series that is summed along with CF2, see 205[:I.J. Thompson and A.R. Barnett, ['Modified Bessel functions I_v and K_v 206 of real order and complex argument to selected accuracy], Computer Physics 207 Communications, vol 47, 245 (1987).] 208 209When ['x] is small (['x] <= 2), CF2 convergence may fail (but CF1 210works very well). The solution here is Temme's series: 211 212[equation mbessel15] 213 214where 215 216[equation mbessel16] 217 218f[sub k] and h[sub k] 219are also computed by recursions (involving gamma functions), but the 220formulas are a little complicated, readers are referred to 221[:N.M. Temme, ['On the numerical evaluation of the modified Bessel function 222 of the third kind], Journal of Computational Physics, vol 19, 324 (1975).] 223Note: Temme's series converge only for |[mu]| <= 1/2. 224 225K[sub [nu]](x) is then calculated from the forward 226recurrence, as is K[sub [nu]+1](x). With these two values and 227f[sub [nu]], the Wronskian yields I[sub [nu]](x) directly. 228 229[endsect] [/section:mbessel Modified Bessel Functions of the First and Second Kinds] 230 231[/ 232 Copyright 2006 John Maddock, Paul A. Bristow and Xiaogang Zhang. 233 Distributed under the Boost Software License, Version 1.0. 234 (See accompanying file LICENSE_1_0.txt or copy at 235 http://www.boost.org/LICENSE_1_0.txt). 236] 237