ceph/src/boost/libs/math/doc/sf/bessel_ik.qbk

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