libs/math/doc/sf_and_dist/bessel_ik.qbk

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