libs/math/doc/sf/bessel_ik.qbk

   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
  22 The functions __cyl_bessel_i and __cyl_bessel_k return the result of the
  23 modified 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
  29 where:
  30
  31 [equation mbessel2]
  32
  33 [equation mbessel3]
  34
  35 The return type of these functions is computed using the __arg_promotion_rules
  36 when T1 and T2 are different types.  The functions are also optimised for the
  37 relatively common case that T1 is an integer.
  38
  39 [optional_policy]
  40
  41 The functions return the result of __domain_error whenever the result is
  42 undefined or complex.  For __cyl_bessel_j this occurs when `x < 0` and v is not
  43 an integer, or when `x == 0` and `v != 0`.  For __cyl_neumann this occurs
  44 when `x <= 0`.
  45
  46 The following graph illustrates the exponential behaviour of I[sub v].
  47
  48 [graph cyl_bessel_i]
  49
  50 The following graph illustrates the exponential decay of K[sub v].
  51
  52 [graph cyl_bessel_k]
  53
  54 [h4 Testing]
  55
  56 There are two sets of test values: spot values calculated using
  57 [@http://functions.wolfram.com functions.wolfram.com],
  58 and a much larger set of tests computed using
  59 a simplified version of this implementation
  60 (with all the special case handling removed).
  61
  62 [h4 Accuracy]
  63
  64 The following tables show how the accuracy of these functions
  65 varies on various platforms, along with comparison to other libraries.
  66 Note that only results for the widest floating-point type on the
  67 system are given, as narrower types have __zero_error.  All values
  68 are relative errors in units of epsilon.  Note that our test suite
  69 includes some fairly extreme inputs which results in most of the worst
  70 problem 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
  80 The following error plot are based on an exhaustive search of the functions domain for I0, I1, K0, and K1,
  81 MSVC-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
 110 The following are handled as special cases first:
 111
 112 When computing I[sub v] for ['x < 0], then [nu] must be an integer
 113 or a domain error occurs.  If [nu] is an integer, then the function is
 114 odd if [nu] is odd and even if [nu] is even, and we can reflect to
 115 ['x > 0].
 116
 117 For I[sub v] with v equal to 0, 1 or 0.5 are handled as special cases.
 118
 119 The 0 and 1 cases use polynomial approximations on
 120 finite and infinite intervals. The approximating forms
 121 are 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"]
 124 by 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
 134 Similarly we have:
 135
 136 [equation bessel22]
 137
 138 [equation bessel23]
 139
 140 [equation bessel24]
 141
 142 [equation bessel25]
 143
 144 The 0.5 case is a simple trigonometric function:
 145
 146 [:I[sub 0.5](x) = sqrt(2 / [pi]x) * sinh(x)]
 147
 148 For K[sub v] with /v/ an integer, the result is calculated using the recurrence relation:
 149
 150 [equation mbessel5]
 151
 152 starting from K[sub 0] and K[sub 1] which are calculated
 153 using rational the approximations above.  These rational approximations are
 154 accurate to around 19 digits, and are therefore only used when T has
 155 no more than 64 binary digits of precision.
 156
 157 When /x/ is small compared to /v/, I[sub v]x is best computed directly from the series:
 158
 159 [equation mbessel17]
 160
 161 In the general case, we first normalize [nu] to \[[^0, [inf]])
 162 with the help of the reflection formulae:
 163
 164 [equation mbessel9]
 165
 166 [equation mbessel10]
 167
 168 Let [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
 170 calculate K[sub [mu]](x) and K[sub [mu]+1](x), and use them to obtain
 171 I[sub [nu]](x) and K[sub [nu]](x).
 172
 173 The 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),]
 176 which needs two continued fractions as well as the Wronskian:
 177
 178 [equation mbessel11]
 179
 180 [equation mbessel12]
 181
 182 [equation mbessel8]
 183
 184 The 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)).]
 187 Their convergence rates depend on ['x], therefore we need
 188 different 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
 194 When ['x] is large (['x] > 2), both continued fractions converge (CF1
 195 may be slow for really large ['x]). K[sub [mu]] and K[sub [mu]+1]
 196 can be calculated by
 197
 198 [equation mbessel13]
 199
 200 where
 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
 209 When ['x] is small (['x] <= 2), CF2 convergence may fail (but CF1
 210 works very well). The solution here is Temme's series:
 211
 212 [equation mbessel15]
 213
 214 where
 215
 216 [equation mbessel16]
 217
 218 f[sub k] and h[sub k]
 219 are also computed by recursions (involving gamma functions), but the
 220 formulas 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).]
 223 Note: Temme's series converge only for |[mu]| <= 1/2.
 224
 225 K[sub [nu]](x) is then calculated from the forward
 226 recurrence, as is K[sub [nu]+1](x). With these two values and
 227 f[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 ]