]>
Commit | Line | Data |
---|---|---|
7c673cae FG |
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 | ] |