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
487
488
489
490
491
492
493
494
495
496
497
498
499
500
501
502
503
504
505
506
507
508
509
510
511
512
513
514
515
516
517
518
519
520
521
522
523
524
525
526
527
528
529
530
531
532
533
534
535
536
537
538
539
540
541
542
543
544
545
546
547
548
549
550
551
552
553
554
555
556
557
558
559
560
561
562
563
564
565
566
567
568
569
570
571
572
573
574
575
576
577
578
579
580
581
582
583
584
585
586
587
588
589
590
591
592
593
594
595
596
597
598
599
600
601
602
603
604
605
606
607
608
609
610
611
612
613
614
615
616
617
618
619
620
621
622
623
624
625
626
627
628
629
630
631
632
633
634
635
636
637
638
639
640
641
642
643
644
645
646
647
648
649
650
651
652
653
654
655
656
657
658
659
660
661
662
663
664
665
666
667
668
669
670
671
672
673
674
675
676
677
678
679
680
681
682
683
684
685
686
687
688
689
690
691
692
693
694
695
696
697
698
|
<html>
<head>
<meta http-equiv="Content-Type" content="text/html; charset=US-ASCII">
<title>Bessel Functions of the First and Second Kinds</title>
<link rel="stylesheet" href="../../math.css" type="text/css">
<meta name="generator" content="DocBook XSL Stylesheets V1.78.1">
<link rel="home" href="../../index.html" title="Math Toolkit 2.1.0">
<link rel="up" href="../bessel.html" title="Bessel Functions">
<link rel="prev" href="bessel_over.html" title="Bessel Function Overview">
<link rel="next" href="bessel_root.html" title="Finding Zeros of 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_over.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="bessel_root.html"><img src="../../../../../../doc/src/images/next.png" alt="Next"></a>
</div>
<div class="section">
<div class="titlepage"><div><div><h3 class="title">
<a name="math_toolkit.bessel.bessel_first"></a><a class="link" href="bessel_first.html" title="Bessel Functions of the First and Second Kinds">Bessel Functions of
the First and Second Kinds</a>
</h3></div></div></div>
<h5>
<a name="math_toolkit.bessel.bessel_first.h0"></a>
<span class="phrase"><a name="math_toolkit.bessel.bessel_first.synopsis"></a></span><a class="link" href="bessel_first.html#math_toolkit.bessel.bessel_first.synopsis">Synopsis</a>
</h5>
<p>
<code class="computeroutput"><span class="preprocessor">#include</span> <span class="special"><</span><span class="identifier">boost</span><span class="special">/</span><span class="identifier">math</span><span class="special">/</span><span class="identifier">special_functions</span><span class="special">/</span><span class="identifier">bessel</span><span class="special">.</span><span class="identifier">hpp</span><span class="special">></span></code>
</p>
<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="../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_j</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="Chapter 14. Policies: Controlling Precision, Error Handling etc">Policy</a><span class="special">></span>
<a class="link" href="../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_j</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="Chapter 14. Policies: Controlling Precision, Error Handling etc">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="../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_neumann</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="Chapter 14. Policies: Controlling Precision, Error Handling etc">Policy</a><span class="special">></span>
<a class="link" href="../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_neumann</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="Chapter 14. Policies: Controlling Precision, Error Handling etc">Policy</a><span class="special">&);</span>
</pre>
<h5>
<a name="math_toolkit.bessel.bessel_first.h1"></a>
<span class="phrase"><a name="math_toolkit.bessel.bessel_first.description"></a></span><a class="link" href="bessel_first.html#math_toolkit.bessel.bessel_first.description">Description</a>
</h5>
<p>
The functions <a class="link" href="bessel_first.html" title="Bessel Functions of the First and Second Kinds">cyl_bessel_j</a>
and <a class="link" href="bessel_first.html" title="Bessel Functions of the First and Second Kinds">cyl_neumann</a> return
the result of the Bessel functions of the first and second kinds respectively:
</p>
<p>
cyl_bessel_j(v, x) = J<sub>v</sub>(x)
</p>
<p>
cyl_neumann(v, x) = Y<sub>v</sub>(x) = N<sub>v</sub>(x)
</p>
<p>
where:
</p>
<p>
<span class="inlinemediaobject"><img src="../../../equations/bessel2.png"></span>
</p>
<p>
<span class="inlinemediaobject"><img src="../../../equations/bessel3.png"></span>
</p>
<p>
The return type of these functions is computed using the <a class="link" href="../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="Chapter 14. Policies: Controlling Precision, Error Handling etc">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="Chapter 14. Policies: Controlling Precision, Error Handling etc">policy
documentation for more details</a>.
</p>
<p>
The functions return the result of <a class="link" href="../error_handling.html#math_toolkit.error_handling.domain_error">domain_error</a>
whenever the result is undefined or complex. For <a class="link" href="bessel_first.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_first.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 cyclic nature of J<sub>v</sub>:
</p>
<p>
<span class="inlinemediaobject"><img src="../../../graphs/cyl_bessel_j.png" align="middle"></span>
</p>
<p>
The following graph shows the behaviour of Y<sub>v</sub>: this is also cyclic for large
<span class="emphasis"><em>x</em></span>, but tends to -∞   for small <span class="emphasis"><em>x</em></span>:
</p>
<p>
<span class="inlinemediaobject"><img src="../../../graphs/cyl_neumann.png" align="middle"></span>
</p>
<h5>
<a name="math_toolkit.bessel.bessel_first.h2"></a>
<span class="phrase"><a name="math_toolkit.bessel.bessel_first.testing"></a></span><a class="link" href="bessel_first.html#math_toolkit.bessel.bessel_first.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.bessel.bessel_first.h3"></a>
<span class="phrase"><a name="math_toolkit.bessel.bessel_first.accuracy"></a></span><a class="link" href="bessel_first.html#math_toolkit.bessel.bessel_first.accuracy">Accuracy</a>
</h5>
<p>
The following tables show how the accuracy of these functions varies on various
platforms, along with comparisons to the <a href="http://www.gnu.org/software/gsl/" target="_top">GSL-1.9</a>
and <a href="http://www.netlib.org/cephes/" target="_top">Cephes</a> libraries.
Note that the cyclic nature of these functions means that they have an infinite
number of irrational roots: in general these functions have arbitrarily large
<span class="emphasis"><em>relative</em></span> errors when the arguments are sufficiently
close to a root. Of course the absolute error in such cases is always small.
Note that only results for the widest floating-point type on the system are
given as narrower types have <a class="link" href="../relative_error.html#math_toolkit.relative_error.zero_error">effectively
zero error</a>. All values are relative errors in units of epsilon.
</p>
<div class="table">
<a name="math_toolkit.bessel.bessel_first.errors_rates_in_cyl_bessel_j"></a><p class="title"><b>Table 6.21. Errors Rates in cyl_bessel_j</b></p>
<div class="table-contents"><table class="table" summary="Errors Rates in cyl_bessel_j">
<colgroup>
<col>
<col>
<col>
<col>
<col>
</colgroup>
<thead><tr>
<th>
<p>
Significand Size
</p>
</th>
<th>
<p>
Platform and Compiler
</p>
</th>
<th>
<p>
J<sub>0</sub>   and J<sub>1</sub>
</p>
</th>
<th>
<p>
J<sub>v</sub>
</p>
</th>
<th>
<p>
J<sub>v</sub>   (large values of x > 1000)
</p>
</th>
</tr></thead>
<tbody>
<tr>
<td>
<p>
53
</p>
</td>
<td>
<p>
Win32 / Visual C++ 8.0
</p>
</td>
<td>
<p>
Peak=2.5 Mean=1.1
</p>
<p>
GSL Peak=6.6
</p>
<p>
<a href="http://www.netlib.org/cephes/" target="_top">Cephes</a> Peak=2.5
Mean=1.1
</p>
</td>
<td>
<p>
Peak=11 Mean=2.2
</p>
<p>
GSL Peak=11
</p>
<p>
<a href="http://www.netlib.org/cephes/" target="_top">Cephes</a> Peak=17
Mean=2.5
</p>
</td>
<td>
<p>
Peak=59 Mean=10
</p>
<p>
GSL Peak=6x10<sup>11</sup>
</p>
<p>
<a href="http://www.netlib.org/cephes/" target="_top">Cephes</a> Peak=2x10<sup>5</sup>
</p>
</td>
</tr>
<tr>
<td>
<p>
64
</p>
</td>
<td>
<p>
Red Hat Linux IA64 / G++ 3.4
</p>
</td>
<td>
<p>
Peak=7 Mean=3
</p>
</td>
<td>
<p>
Peak=117 Mean=10
</p>
</td>
<td>
<p>
Peak=2x10<sup>4</sup>   Mean=6x10<sup>3</sup>
</p>
</td>
</tr>
<tr>
<td>
<p>
64
</p>
</td>
<td>
<p>
SUSE Linux AMD64 / G++ 4.1
</p>
</td>
<td>
<p>
Peak=7 Mean=3
</p>
</td>
<td>
<p>
Peak=400 Mean=40
</p>
</td>
<td>
<p>
Peak=2x10<sup>4</sup>   Mean=1x10<sup>4</sup>
</p>
</td>
</tr>
<tr>
<td>
<p>
113
</p>
</td>
<td>
<p>
HP-UX / HP aCC 6
</p>
</td>
<td>
<p>
Peak=14 Mean=6
</p>
</td>
<td>
<p>
Peak=29 Mean=3
</p>
</td>
<td>
<p>
Peak=2700 Mean=450
</p>
</td>
</tr>
</tbody>
</table></div>
</div>
<br class="table-break"><div class="table">
<a name="math_toolkit.bessel.bessel_first.errors_rates_in_cyl_neumann"></a><p class="title"><b>Table 6.22. Errors Rates in cyl_neumann</b></p>
<div class="table-contents"><table class="table" summary="Errors Rates in cyl_neumann">
<colgroup>
<col>
<col>
<col>
<col>
<col>
</colgroup>
<thead><tr>
<th>
<p>
Significand Size
</p>
</th>
<th>
<p>
Platform and Compiler
</p>
</th>
<th>
<p>
Y<sub>0</sub>   and Y<sub>1</sub>
</p>
</th>
<th>
<p>
Y<sub>n</sub> (integer orders)
</p>
</th>
<th>
<p>
Y<sub>v</sub> (fractional orders)
</p>
</th>
</tr></thead>
<tbody>
<tr>
<td>
<p>
53
</p>
</td>
<td>
<p>
Win32 / Visual C++ 8.0
</p>
</td>
<td>
<p>
Peak=4.7 Mean=1.7
</p>
<p>
GSL Peak=34 Mean=9
</p>
<p>
<a href="http://www.netlib.org/cephes/" target="_top">Cephes</a> Peak=330
Mean=54
</p>
</td>
<td>
<p>
Peak=117 Mean=10
</p>
<p>
GSL Peak=500 Mean=54
</p>
<p>
<a href="http://www.netlib.org/cephes/" target="_top">Cephes</a> Peak=923
Mean=83
</p>
</td>
<td>
<p>
Peak=800 Mean=40
</p>
<p>
GSL Peak=1.4x10<sup>6</sup>   Mean=7x10<sup>4</sup>  
</p>
<p>
<a href="http://www.netlib.org/cephes/" target="_top">Cephes</a> Peak=+INF
</p>
</td>
</tr>
<tr>
<td>
<p>
64
</p>
</td>
<td>
<p>
Red Hat Linux IA64 / G++ 3.4
</p>
</td>
<td>
<p>
Peak=470 Mean=56
</p>
</td>
<td>
<p>
Peak=843 Mean=51
</p>
</td>
<td>
<p>
Peak=741 Mean=51
</p>
</td>
</tr>
<tr>
<td>
<p>
64
</p>
</td>
<td>
<p>
SUSE Linux AMD64 / G++ 4.1
</p>
</td>
<td>
<p>
Peak=1300 Mean=424
</p>
</td>
<td>
<p>
Peak=2x10<sup>4</sup>   Mean=8x10<sup>3</sup>
</p>
</td>
<td>
<p>
Peak=1x10<sup>5</sup>   Mean=6x10<sup>3</sup>
</p>
</td>
</tr>
<tr>
<td>
<p>
113
</p>
</td>
<td>
<p>
HP-UX / HP aCC 6
</p>
</td>
<td>
<p>
Peak=180 Mean=63
</p>
</td>
<td>
<p>
Peak=340 Mean=150
</p>
</td>
<td>
<p>
Peak=2x10<sup>4</sup>   Mean=1200
</p>
</td>
</tr>
</tbody>
</table></div>
</div>
<br class="table-break"><p>
Note that for large <span class="emphasis"><em>x</em></span> these functions are largely dependent
on the accuracy of the <code class="computeroutput"><span class="identifier">std</span><span class="special">::</span><span class="identifier">sin</span></code> and
<code class="computeroutput"><span class="identifier">std</span><span class="special">::</span><span class="identifier">cos</span></code> functions.
</p>
<p>
Comparison to GSL and <a href="http://www.netlib.org/cephes/" target="_top">Cephes</a>
is interesting: both <a href="http://www.netlib.org/cephes/" target="_top">Cephes</a>
and this library optimise the integer order case - leading to identical results
- simply using the general case is for the most part slightly more accurate
though, as noted by the better accuracy of GSL in the integer argument cases.
This implementation tends to perform much better when the arguments become
large, <a href="http://www.netlib.org/cephes/" target="_top">Cephes</a> in particular
produces some remarkably inaccurate results with some of the test data (no
significant figures correct), and even GSL performs badly with some inputs
to J<sub>v</sub>. Note that by way of double-checking these results, the worst performing
<a href="http://www.netlib.org/cephes/" target="_top">Cephes</a> and GSL cases were
recomputed using <a href="http://functions.wolfram.com" target="_top">functions.wolfram.com</a>,
and the result checked against our test data: no errors in the test data
were found.
</p>
<h5>
<a name="math_toolkit.bessel.bessel_first.h4"></a>
<span class="phrase"><a name="math_toolkit.bessel.bessel_first.implementation"></a></span><a class="link" href="bessel_first.html#math_toolkit.bessel.bessel_first.implementation">Implementation</a>
</h5>
<p>
The implementation is mostly about filtering off various special cases:
</p>
<p>
When <span class="emphasis"><em>x</em></span> is negative, then the order <span class="emphasis"><em>v</em></span>
must be an integer or the result is a domain error. If the order is an integer
then the function is odd for odd orders and even for even orders, so we reflect
to <span class="emphasis"><em>x > 0</em></span>.
</p>
<p>
When the order <span class="emphasis"><em>v</em></span> is negative then the reflection formulae
can be used to move to <span class="emphasis"><em>v > 0</em></span>:
</p>
<p>
<span class="inlinemediaobject"><img src="../../../equations/bessel9.png"></span>
</p>
<p>
<span class="inlinemediaobject"><img src="../../../equations/bessel10.png"></span>
</p>
<p>
Note that if the order is an integer, then these formulae reduce to:
</p>
<p>
J<sub>-n</sub> = (-1)<sup>n</sup>J<sub>n</sub>
</p>
<p>
Y<sub>-n</sub> = (-1)<sup>n</sup>Y<sub>n</sub>
</p>
<p>
However, in general, a negative order implies that we will need to compute
both J and Y.
</p>
<p>
When <span class="emphasis"><em>x</em></span> is large compared to the order <span class="emphasis"><em>v</em></span>
then the asymptotic expansions for large <span class="emphasis"><em>x</em></span> in M. Abramowitz
and I.A. Stegun, <span class="emphasis"><em>Handbook of Mathematical Functions</em></span>
9.2.19 are used (these were found to be more reliable than those in A&S
9.2.5).
</p>
<p>
When the order <span class="emphasis"><em>v</em></span> is an integer the method first relates
the result to J<sub>0</sub>, J<sub>1</sub>, Y<sub>0</sub>   and Y<sub>1</sub>   using either forwards or backwards recurrence
(Miller's algorithm) depending upon which is stable. The values for J<sub>0</sub>, J<sub>1</sub>,
Y<sub>0</sub>   and Y<sub>1</sub>   are calculated using the rational minimax approximations on root-bracketing
intervals for small <span class="emphasis"><em>|x|</em></span> and Hankel asymptotic expansion
for large <span class="emphasis"><em>|x|</em></span>. The coefficients are from:
</p>
<p>
W.J. Cody, <span class="emphasis"><em>ALGORITHM 715: SPECFUN - A Portable FORTRAN Package
of Special Function Routines and Test Drivers</em></span>, ACM Transactions
on Mathematical Software, vol 19, 22 (1993).
</p>
<p>
and
</p>
<p>
J.F. Hart et al, <span class="emphasis"><em>Computer Approximations</em></span>, John Wiley
& Sons, New York, 1968.
</p>
<p>
These approximations are accurate to around 19 decimal digits: therefore
these methods are not used when type T has more than 64 binary digits.
</p>
<p>
When <span class="emphasis"><em>x</em></span> is smaller than machine epsilon then the following
approximations for Y<sub>0</sub>(x), Y<sub>1</sub>(x), Y<sub>2</sub>(x) and Y<sub>n</sub>(x) can be used (see: <a href="http://functions.wolfram.com/03.03.06.0037.01" target="_top">http://functions.wolfram.com/03.03.06.0037.01</a>,
<a href="http://functions.wolfram.com/03.03.06.0038.01" target="_top">http://functions.wolfram.com/03.03.06.0038.01</a>,
<a href="http://functions.wolfram.com/03.03.06.0039.01" target="_top">http://functions.wolfram.com/03.03.06.0039.01</a>
and <a href="http://functions.wolfram.com/03.03.06.0040.01" target="_top">http://functions.wolfram.com/03.03.06.0040.01</a>):
</p>
<p>
<span class="inlinemediaobject"><img src="../../../equations/bessel_y0_small_z.png"></span>
</p>
<p>
<span class="inlinemediaobject"><img src="../../../equations/bessel_y1_small_z.png"></span>
</p>
<p>
<span class="inlinemediaobject"><img src="../../../equations/bessel_y2_small_z.png"></span>
</p>
<p>
<span class="inlinemediaobject"><img src="../../../equations/bessel_yn_small_z.png"></span>
</p>
<p>
When <span class="emphasis"><em>x</em></span> is small compared to <span class="emphasis"><em>v</em></span> and
<span class="emphasis"><em>v</em></span> is not an integer, then the following series approximation
can be used for Y<sub>v</sub>(x), this is also an area where other approximations are
often too slow to converge to be used (see <a href="http://functions.wolfram.com/03.03.06.0034.01" target="_top">http://functions.wolfram.com/03.03.06.0034.01</a>):
</p>
<p>
<span class="inlinemediaobject"><img src="../../../equations/bessel_yv_small_z.png"></span>
</p>
<p>
When <span class="emphasis"><em>x</em></span> is small compared to <span class="emphasis"><em>v</em></span>,
J<sub>v</sub>x   is best computed directly from the series:
</p>
<p>
<span class="inlinemediaobject"><img src="../../../equations/bessel2.png"></span>
</p>
<p>
In the general case we compute J<sub>v</sub>   and Y<sub>v</sub>   simultaneously.
</p>
<p>
To get the initial values, 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 J<sub>μ</sub>(x), J<sub>μ+1</sub>(x), Y<sub>μ</sub>(x), Y<sub>μ+1</sub>(x) and use them to obtain J<sub>ν</sub>(x), Y<sub>ν</sub>(x).
</p>
<p>
The algorithm is called Steed's method, which needs two continued fractions
as well as the Wronskian:
</p>
<p>
<span class="inlinemediaobject"><img src="../../../equations/bessel8.png"></span>
</p>
<p>
<span class="inlinemediaobject"><img src="../../../equations/bessel11.png"></span>
</p>
<p>
<span class="inlinemediaobject"><img src="../../../equations/bessel12.png"></span>
</p>
<p>
See: F.S. Acton, <span class="emphasis"><em>Numerical Methods that Work</em></span>, The Mathematical
Association of America, Washington, 1997.
</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>).
J<sub>μ</sub>, J<sub>μ+1</sub>, Y<sub>μ</sub>, Y<sub>μ+1</sub> can be calculated by
</p>
<p>
<span class="inlinemediaobject"><img src="../../../equations/bessel13.png"></span>
</p>
<p>
where
</p>
<p>
<span class="inlinemediaobject"><img src="../../../equations/bessel14.png"></span>
</p>
<p>
J<sub>ν</sub> and Y<sub>μ</sub> are then calculated using backward (Miller's algorithm) and forward
recurrence respectively.
</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/bessel15.png"></span>
</p>
<p>
where
</p>
<p>
<span class="inlinemediaobject"><img src="../../../equations/bessel16.png"></span>
</p>
<p>
g<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 refered to N.M. Temme,
<span class="emphasis"><em>On the numerical evaluation of the ordinary Bessel function of
the second kind</em></span>, Journal of Computational Physics, vol 21, 343
(1976). Note Temme's series converge only for |μ| <= 1/2.
</p>
<p>
As the previous case, Y<sub>ν</sub>   is calculated from the forward recurrence, so is Y<sub>ν+1</sub>.
With these two values and f<sub>ν</sub>, the Wronskian yields J<sub>ν</sub>(x) directly without backward
recurrence.
</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, 2012-2014 Nikhar Agrawal,
Anton Bikineev, Paul A. Bristow, Marco Guazzone, Christopher Kormanyos, Hubert
Holin, Bruno Lalande, John Maddock, Johan Råde, Gautam Sewani, Benjamin Sobotta,
Thijs van den Berg, Daryle Walker and Xiaogang Zhang<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_over.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="bessel_root.html"><img src="../../../../../../doc/src/images/next.png" alt="Next"></a>
</div>
</body>
</html>
|
