1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
419
420
421
422
423
424
425
426
427
428
429
430
431
432
433
434
435
436
437
438
439
440
441
442
443
444
445
446
447
448
449
450
451
452
453
454
455
456
457
458
459
460
461
462
463
464
465
466
467
468
469
470
471
472
473
474
475
476
477
478
479
480
481
482
483
484
485
486
|
<html>
<head>
<meta http-equiv="Content-Type" content="text/html; charset=US-ASCII">
<title>Modified Bessel Functions of the First and Second Kinds</title>
<link rel="stylesheet" href="../../../../../../../../doc/src/boostbook.css" type="text/css">
<meta name="generator" content="DocBook XSL Stylesheets V1.76.1">
<link rel="home" href="../../../index.html" title="Math Toolkit">
<link rel="up" href="../bessel.html" title="Bessel Functions">
<link rel="prev" href="bessel.html" title="Bessel Functions of the First and Second Kinds">
<link rel="next" href="sph_bessel.html" title="Spherical Bessel Functions of the First and Second Kinds">
</head>
<body bgcolor="white" text="black" link="#0000FF" vlink="#840084" alink="#0000FF">
<table cellpadding="2" width="100%"><tr>
<td valign="top"><img alt="Boost C++ Libraries" width="277" height="86" src="../../../../../../../../boost.png"></td>
<td align="center"><a href="../../../../../../../../index.html">Home</a></td>
<td align="center"><a href="../../../../../../../../libs/libraries.htm">Libraries</a></td>
<td align="center"><a href="http://www.boost.org/users/people.html">People</a></td>
<td align="center"><a href="http://www.boost.org/users/faq.html">FAQ</a></td>
<td align="center"><a href="../../../../../../../../more/index.htm">More</a></td>
</tr></table>
<hr>
<div class="spirit-nav">
<a accesskey="p" href="bessel.html"><img src="../../../../../../../../doc/src/images/prev.png" alt="Prev"></a><a accesskey="u" href="../bessel.html"><img src="../../../../../../../../doc/src/images/up.png" alt="Up"></a><a accesskey="h" href="../../../index.html"><img src="../../../../../../../../doc/src/images/home.png" alt="Home"></a><a accesskey="n" href="sph_bessel.html"><img src="../../../../../../../../doc/src/images/next.png" alt="Next"></a>
</div>
<div class="section math_toolkit_special_bessel_mbessel">
<div class="titlepage"><div><div><h4 class="title">
<a name="math_toolkit.special.bessel.mbessel"></a><a class="link" href="mbessel.html" title="Modified Bessel Functions of the First and Second Kinds">Modified Bessel
Functions of the First and Second Kinds</a>
</h4></div></div></div>
<h5>
<a name="math_toolkit.special.bessel.mbessel.h0"></a>
<span><a name="math_toolkit.special.bessel.mbessel.synopsis"></a></span><a class="link" href="mbessel.html#math_toolkit.special.bessel.mbessel.synopsis">Synopsis</a>
</h5>
<pre class="programlisting"><span class="keyword">template</span> <span class="special"><</span><span class="keyword">class</span> <span class="identifier">T1</span><span class="special">,</span> <span class="keyword">class</span> <span class="identifier">T2</span><span class="special">></span>
<a class="link" href="../../main_overview/result_type.html" title="Calculation of the Type of the Result"><span class="emphasis"><em>calculated-result-type</em></span></a> <span class="identifier">cyl_bessel_i</span><span class="special">(</span><span class="identifier">T1</span> <span class="identifier">v</span><span class="special">,</span> <span class="identifier">T2</span> <span class="identifier">x</span><span class="special">);</span>
<span class="keyword">template</span> <span class="special"><</span><span class="keyword">class</span> <span class="identifier">T1</span><span class="special">,</span> <span class="keyword">class</span> <span class="identifier">T2</span><span class="special">,</span> <span class="keyword">class</span> <a class="link" href="../../policy.html" title="Policies">Policy</a><span class="special">></span>
<a class="link" href="../../main_overview/result_type.html" title="Calculation of the Type of the Result"><span class="emphasis"><em>calculated-result-type</em></span></a> <span class="identifier">cyl_bessel_i</span><span class="special">(</span><span class="identifier">T1</span> <span class="identifier">v</span><span class="special">,</span> <span class="identifier">T2</span> <span class="identifier">x</span><span class="special">,</span> <span class="keyword">const</span> <a class="link" href="../../policy.html" title="Policies">Policy</a><span class="special">&);</span>
<span class="keyword">template</span> <span class="special"><</span><span class="keyword">class</span> <span class="identifier">T1</span><span class="special">,</span> <span class="keyword">class</span> <span class="identifier">T2</span><span class="special">></span>
<a class="link" href="../../main_overview/result_type.html" title="Calculation of the Type of the Result"><span class="emphasis"><em>calculated-result-type</em></span></a> <span class="identifier">cyl_bessel_k</span><span class="special">(</span><span class="identifier">T1</span> <span class="identifier">v</span><span class="special">,</span> <span class="identifier">T2</span> <span class="identifier">x</span><span class="special">);</span>
<span class="keyword">template</span> <span class="special"><</span><span class="keyword">class</span> <span class="identifier">T1</span><span class="special">,</span> <span class="keyword">class</span> <span class="identifier">T2</span><span class="special">,</span> <span class="keyword">class</span> <a class="link" href="../../policy.html" title="Policies">Policy</a><span class="special">></span>
<a class="link" href="../../main_overview/result_type.html" title="Calculation of the Type of the Result"><span class="emphasis"><em>calculated-result-type</em></span></a> <span class="identifier">cyl_bessel_k</span><span class="special">(</span><span class="identifier">T1</span> <span class="identifier">v</span><span class="special">,</span> <span class="identifier">T2</span> <span class="identifier">x</span><span class="special">,</span> <span class="keyword">const</span> <a class="link" href="../../policy.html" title="Policies">Policy</a><span class="special">&);</span>
</pre>
<h5>
<a name="math_toolkit.special.bessel.mbessel.h1"></a>
<span><a name="math_toolkit.special.bessel.mbessel.description"></a></span><a class="link" href="mbessel.html#math_toolkit.special.bessel.mbessel.description">Description</a>
</h5>
<p>
The functions <a class="link" href="mbessel.html" title="Modified Bessel Functions of the First and Second Kinds">cyl_bessel_i</a>
and <a class="link" href="mbessel.html" title="Modified Bessel Functions of the First and Second Kinds">cyl_bessel_k</a>
return the result of the modified Bessel functions of the first and second
kind respectively:
</p>
<p>
cyl_bessel_i(v, x) = I<sub>v</sub>(x)
</p>
<p>
cyl_bessel_k(v, x) = K<sub>v</sub>(x)
</p>
<p>
where:
</p>
<p>
<span class="inlinemediaobject"><img src="../../../../equations/mbessel2.png"></span>
</p>
<p>
<span class="inlinemediaobject"><img src="../../../../equations/mbessel3.png"></span>
</p>
<p>
The return type of these functions is computed using the <a class="link" href="../../main_overview/result_type.html" title="Calculation of the Type of the Result"><span class="emphasis"><em>result
type calculation rules</em></span></a> when T1 and T2 are different types.
The functions are also optimised for the relatively common case that T1
is an integer.
</p>
<p>
The final <a class="link" href="../../policy.html" title="Policies">Policy</a> argument is
optional and can be used to control the behaviour of the function: how
it handles errors, what level of precision to use etc. Refer to the <a class="link" href="../../policy.html" title="Policies">policy documentation for more details</a>.
</p>
<p>
The functions return the result of <a class="link" href="../../main_overview/error_handling.html#domain_error">domain_error</a>
whenever the result is undefined or complex. For <a class="link" href="bessel.html" title="Bessel Functions of the First and Second Kinds">cyl_bessel_j</a>
this occurs when <code class="computeroutput"><span class="identifier">x</span> <span class="special"><</span>
<span class="number">0</span></code> and v is not an integer, or when
<code class="computeroutput"><span class="identifier">x</span> <span class="special">==</span>
<span class="number">0</span></code> and <code class="computeroutput"><span class="identifier">v</span>
<span class="special">!=</span> <span class="number">0</span></code>.
For <a class="link" href="bessel.html" title="Bessel Functions of the First and Second Kinds">cyl_neumann</a>
this occurs when <code class="computeroutput"><span class="identifier">x</span> <span class="special"><=</span>
<span class="number">0</span></code>.
</p>
<p>
The following graph illustrates the exponential behaviour of I<sub>v</sub>.
</p>
<p>
<span class="inlinemediaobject"><img src="../../../../graphs/cyl_bessel_i.png" align="middle"></span>
</p>
<p>
The following graph illustrates the exponential decay of K<sub>v</sub>.
</p>
<p>
<span class="inlinemediaobject"><img src="../../../../graphs/cyl_bessel_k.png" align="middle"></span>
</p>
<h5>
<a name="math_toolkit.special.bessel.mbessel.h2"></a>
<span><a name="math_toolkit.special.bessel.mbessel.testing"></a></span><a class="link" href="mbessel.html#math_toolkit.special.bessel.mbessel.testing">Testing</a>
</h5>
<p>
There are two sets of test values: spot values calculated using <a href="http://functions.wolfram.com" target="_top">functions.wolfram.com</a>, and a
much larger set of tests computed using a simplified version of this implementation
(with all the special case handling removed).
</p>
<h5>
<a name="math_toolkit.special.bessel.mbessel.h3"></a>
<span><a name="math_toolkit.special.bessel.mbessel.accuracy"></a></span><a class="link" href="mbessel.html#math_toolkit.special.bessel.mbessel.accuracy">Accuracy</a>
</h5>
<p>
The following tables show how the accuracy of these functions varies on
various platforms, along with a comparison to the <a href="http://www.gnu.org/software/gsl/" target="_top">GSL-1.9</a>
library. Note that only results for the widest floating-point type on the
system are given, as narrower types have <a class="link" href="../../backgrounders/relative_error.html#zero_error">effectively
zero error</a>. All values are relative errors in units of epsilon.
</p>
<div class="table">
<a name="math_toolkit.special.bessel.mbessel.errors_rates_in_cyl_bessel_i"></a><p class="title"><b>Table 41. Errors Rates in cyl_bessel_i</b></p>
<div class="table-contents"><table class="table" summary="Errors Rates in cyl_bessel_i">
<colgroup>
<col>
<col>
<col>
</colgroup>
<thead><tr>
<th>
<p>
Significand Size
</p>
</th>
<th>
<p>
Platform and Compiler
</p>
</th>
<th>
<p>
I<sub>v</sub>
</p>
</th>
</tr></thead>
<tbody>
<tr>
<td>
<p>
53
</p>
</td>
<td>
<p>
Win32 / Visual C++ 8.0
</p>
</td>
<td>
<p>
Peak=10 Mean=3.4 GSL Peak=6000
</p>
</td>
</tr>
<tr>
<td>
<p>
64
</p>
</td>
<td>
<p>
Red Hat Linux IA64 / G++ 3.4
</p>
</td>
<td>
<p>
Peak=11 Mean=3
</p>
</td>
</tr>
<tr>
<td>
<p>
64
</p>
</td>
<td>
<p>
SUSE Linux AMD64 / G++ 4.1
</p>
</td>
<td>
<p>
Peak=11 Mean=4
</p>
</td>
</tr>
<tr>
<td>
<p>
113
</p>
</td>
<td>
<p>
HP-UX / HP aCC 6
</p>
</td>
<td>
<p>
Peak=15 Mean=4
</p>
</td>
</tr>
</tbody>
</table></div>
</div>
<br class="table-break"><div class="table">
<a name="math_toolkit.special.bessel.mbessel.errors_rates_in_cyl_bessel_k"></a><p class="title"><b>Table 42. Errors Rates in cyl_bessel_k</b></p>
<div class="table-contents"><table class="table" summary="Errors Rates in cyl_bessel_k">
<colgroup>
<col>
<col>
<col>
</colgroup>
<thead><tr>
<th>
<p>
Significand Size
</p>
</th>
<th>
<p>
Platform and Compiler
</p>
</th>
<th>
<p>
K<sub>v</sub>
</p>
</th>
</tr></thead>
<tbody>
<tr>
<td>
<p>
53
</p>
</td>
<td>
<p>
Win32 / Visual C++ 8.0
</p>
</td>
<td>
<p>
Peak=9 Mean=2
</p>
<p>
GSL Peak=9
</p>
</td>
</tr>
<tr>
<td>
<p>
64
</p>
</td>
<td>
<p>
Red Hat Linux IA64 / G++ 3.4
</p>
</td>
<td>
<p>
Peak=10 Mean=2
</p>
</td>
</tr>
<tr>
<td>
<p>
64
</p>
</td>
<td>
<p>
SUSE Linux AMD64 / G++ 4.1
</p>
</td>
<td>
<p>
Peak=10 Mean=2
</p>
</td>
</tr>
<tr>
<td>
<p>
113
</p>
</td>
<td>
<p>
HP-UX / HP aCC 6
</p>
</td>
<td>
<p>
Peak=12 Mean=5
</p>
</td>
</tr>
</tbody>
</table></div>
</div>
<br class="table-break"><h5>
<a name="math_toolkit.special.bessel.mbessel.h4"></a>
<span><a name="math_toolkit.special.bessel.mbessel.implementation"></a></span><a class="link" href="mbessel.html#math_toolkit.special.bessel.mbessel.implementation">Implementation</a>
</h5>
<p>
The following are handled as special cases first:
</p>
<p>
When computing I<sub>v</sub>   for <span class="emphasis"><em>x < 0</em></span>, then ν   must be an integer
or a domain error occurs. If ν   is an integer, then the function is odd if
ν   is odd and even if ν   is even, and we can reflect to <span class="emphasis"><em>x > 0</em></span>.
</p>
<p>
For I<sub>v</sub>   with v equal to 0, 1 or 0.5 are handled as special cases.
</p>
<p>
The 0 and 1 cases use minimax rational approximations on finite and infinite
intervals. The coefficients are from:
</p>
<div class="itemizedlist"><ul class="itemizedlist" type="disc">
<li class="listitem">
J.M. Blair and C.A. Edwards, <span class="emphasis"><em>Stable rational minimax approximations
to the modified Bessel functions I_0(x) and I_1(x)</em></span>, Atomic
Energy of Canada Limited Report 4928, Chalk River, 1974.
</li>
<li class="listitem">
S. Moshier, <span class="emphasis"><em>Methods and Programs for Mathematical Functions</em></span>,
Ellis Horwood Ltd, Chichester, 1989.
</li>
</ul></div>
<p>
While the 0.5 case is a simple trigonometric function:
</p>
<p>
I<sub>0.5</sub>(x) = sqrt(2 / πx) * sinh(x)
</p>
<p>
For K<sub>v</sub>   with <span class="emphasis"><em>v</em></span> an integer, the result is calculated
using the recurrence relation:
</p>
<p>
<span class="inlinemediaobject"><img src="../../../../equations/mbessel5.png"></span>
</p>
<p>
starting from K<sub>0</sub>   and K<sub>1</sub>   which are calculated using rational the approximations
above. These rational approximations are accurate to around 19 digits,
and are therefore only used when T has no more than 64 binary digits of
precision.
</p>
<p>
When <span class="emphasis"><em>x</em></span> is small compared to <span class="emphasis"><em>v</em></span>,
I<sub>v</sub>x   is best computed directly from the series:
</p>
<p>
<span class="inlinemediaobject"><img src="../../../../equations/mbessel17.png"></span>
</p>
<p>
In the general case, we first normalize ν   to [<code class="literal">0, [inf</code>])
with the help of the reflection formulae:
</p>
<p>
<span class="inlinemediaobject"><img src="../../../../equations/mbessel9.png"></span>
</p>
<p>
<span class="inlinemediaobject"><img src="../../../../equations/mbessel10.png"></span>
</p>
<p>
Let μ   = ν - floor(ν + 1/2), then μ   is the fractional part of ν   such that |μ| <=
1/2 (we need this for convergence later). The idea is to calculate K<sub>μ</sub>(x)
and K<sub>μ+1</sub>(x), and use them to obtain I<sub>ν</sub>(x) and K<sub>ν</sub>(x).
</p>
<p>
The algorithm is proposed by Temme in N.M. Temme, <span class="emphasis"><em>On the numerical
evaluation of the modified bessel function of the third kind</em></span>,
Journal of Computational Physics, vol 19, 324 (1975), which needs two continued
fractions as well as the Wronskian:
</p>
<p>
<span class="inlinemediaobject"><img src="../../../../equations/mbessel11.png"></span>
</p>
<p>
<span class="inlinemediaobject"><img src="../../../../equations/mbessel12.png"></span>
</p>
<p>
<span class="inlinemediaobject"><img src="../../../../equations/mbessel8.png"></span>
</p>
<p>
The continued fractions are computed using the modified Lentz's method
(W.J. Lentz, <span class="emphasis"><em>Generating Bessel functions in Mie scattering calculations
using continued fractions</em></span>, Applied Optics, vol 15, 668 (1976)).
Their convergence rates depend on <span class="emphasis"><em>x</em></span>, therefore we
need different strategies for large <span class="emphasis"><em>x</em></span> and small <span class="emphasis"><em>x</em></span>.
</p>
<p>
<span class="emphasis"><em>x > v</em></span>, CF1 needs O(<span class="emphasis"><em>x</em></span>) iterations
to converge, CF2 converges rapidly.
</p>
<p>
<span class="emphasis"><em>x <= v</em></span>, CF1 converges rapidly, CF2 fails to converge
when <span class="emphasis"><em>x</em></span> <code class="literal">-></code> 0.
</p>
<p>
When <span class="emphasis"><em>x</em></span> is large (<span class="emphasis"><em>x</em></span> > 2), both
continued fractions converge (CF1 may be slow for really large <span class="emphasis"><em>x</em></span>).
K<sub>μ</sub>   and K<sub>μ+1</sub>  
can be calculated by
</p>
<p>
<span class="inlinemediaobject"><img src="../../../../equations/mbessel13.png"></span>
</p>
<p>
where
</p>
<p>
<span class="inlinemediaobject"><img src="../../../../equations/mbessel14.png"></span>
</p>
<p>
<span class="emphasis"><em>S</em></span> is also a series that is summed along with CF2,
see I.J. Thompson and A.R. Barnett, <span class="emphasis"><em>Modified Bessel functions
I_v and K_v of real order and complex argument to selected accuracy</em></span>,
Computer Physics Communications, vol 47, 245 (1987).
</p>
<p>
When <span class="emphasis"><em>x</em></span> is small (<span class="emphasis"><em>x</em></span> <= 2),
CF2 convergence may fail (but CF1 works very well). The solution here is
Temme's series:
</p>
<p>
<span class="inlinemediaobject"><img src="../../../../equations/mbessel15.png"></span>
</p>
<p>
where
</p>
<p>
<span class="inlinemediaobject"><img src="../../../../equations/mbessel16.png"></span>
</p>
<p>
f<sub>k</sub>   and h<sub>k</sub>  
are also computed by recursions (involving gamma functions), but
the formulas are a little complicated, readers are referred to N.M. Temme,
<span class="emphasis"><em>On the numerical evaluation of the modified Bessel function of
the third kind</em></span>, Journal of Computational Physics, vol 19, 324
(1975). Note: Temme's series converge only for |μ| <= 1/2.
</p>
<p>
K<sub>ν</sub>(x) is then calculated from the forward recurrence, as is K<sub>ν+1</sub>(x). With
these two values and f<sub>ν</sub>, the Wronskian yields I<sub>ν</sub>(x) directly.
</p>
</div>
<table xmlns:rev="http://www.cs.rpi.edu/~gregod/boost/tools/doc/revision" width="100%"><tr>
<td align="left"></td>
<td align="right"><div class="copyright-footer">Copyright © 2006-2010 John Maddock, Paul A. Bristow, Hubert Holin, Xiaogang Zhang, Bruno
Lalande, Johan Råde, Gautam Sewani, Thijs van den Berg and Benjamin Sobotta<p>
Distributed under the Boost Software License, Version 1.0. (See accompanying
file LICENSE_1_0.txt or copy at <a href="http://www.boost.org/LICENSE_1_0.txt" target="_top">http://www.boost.org/LICENSE_1_0.txt</a>)
</p>
</div></td>
</tr></table>
<hr>
<div class="spirit-nav">
<a accesskey="p" href="bessel.html"><img src="../../../../../../../../doc/src/images/prev.png" alt="Prev"></a><a accesskey="u" href="../bessel.html"><img src="../../../../../../../../doc/src/images/up.png" alt="Up"></a><a accesskey="h" href="../../../index.html"><img src="../../../../../../../../doc/src/images/home.png" alt="Home"></a><a accesskey="n" href="sph_bessel.html"><img src="../../../../../../../../doc/src/images/next.png" alt="Next"></a>
</div>
</body>
</html>
|
