2 [section:bessel_first Bessel Functions of the First and Second Kinds]
6 `#include <boost/math/special_functions/bessel.hpp>`
8 template <class T1, class T2>
9 ``__sf_result`` cyl_bessel_j(T1 v, T2 x);
11 template <class T1, class T2, class ``__Policy``>
12 ``__sf_result`` cyl_bessel_j(T1 v, T2 x, const ``__Policy``&);
14 template <class T1, class T2>
15 ``__sf_result`` cyl_neumann(T1 v, T2 x);
17 template <class T1, class T2, class ``__Policy``>
18 ``__sf_result`` cyl_neumann(T1 v, T2 x, const ``__Policy``&);
23 The functions __cyl_bessel_j and __cyl_neumann return the result of the
24 Bessel functions of the first and second kinds respectively:
26 cyl_bessel_j(v, x) = J[sub v](x)
28 cyl_neumann(v, x) = Y[sub v](x) = N[sub v](x)
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.
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
47 The following graph illustrates the cyclic nature of J[sub v]:
51 The following graph shows the behaviour of Y[sub v]: this is also
52 cyclic for large /x/, but tends to -[infin][space] for small /x/:
58 There are two sets of test values: spot values calculated using
59 [@http://functions.wolfram.com functions.wolfram.com],
60 and a much larger set of tests computed using
61 a simplified version of this implementation
62 (with all the special case handling removed).
66 The following tables show how the accuracy of these functions
67 varies on various platforms, along with comparisons to other
68 libraries. Note that the cyclic nature of these
69 functions means that they have an infinite number of irrational
70 roots: in general these functions have arbitrarily large /relative/
71 errors when the arguments are sufficiently close to a root. Of
72 course the absolute error in such cases is always small.
73 Note that only results for the widest floating-point type on the
74 system are given as narrower types have __zero_error. All values
75 are relative errors in units of epsilon. Most of the gross errors
76 exhibited by other libraries occur for very large arguments - you will
77 need to drill down into the actual program output if you need more
80 [table_cyl_bessel_j_integer_orders_]
84 [table_cyl_neumann_integer_orders_]
88 Note that for large /x/ these functions are largely dependent on
89 the accuracy of the `std::sin` and `std::cos` functions.
91 Comparison to GSL and __cephes is interesting: both __cephes and this library optimise
92 the integer order case - leading to identical results - simply using the general
93 case is for the most part slightly more accurate though, as noted by the
94 better accuracy of GSL in the integer argument cases. This implementation tends to
95 perform much better when the arguments become large, __cephes in particular produces
96 some remarkably inaccurate results with some of the test data (no significant figures
97 correct), and even GSL performs badly with some inputs to J[sub v]. Note that
98 by way of double-checking these results, the worst performing __cephes and GSL cases
99 were recomputed using [@http://functions.wolfram.com functions.wolfram.com],
100 and the result checked against our test data: no errors in the test data were found.
104 The implementation is mostly about filtering off various special cases:
106 When /x/ is negative, then the order /v/ must be an integer or the
107 result is a domain error. If the order is an integer then the function
108 is odd for odd orders and even for even orders, so we reflect to /x > 0/.
110 When the order /v/ is negative then the reflection formulae can be used to
117 Note that if the order is an integer, then these formulae reduce to:
119 J[sub -n] = (-1)[super n]J[sub n]
121 Y[sub -n] = (-1)[super n]Y[sub n]
123 However, in general, a negative order implies that we will need to compute
126 When /x/ is large compared to the order /v/ then the asymptotic expansions
127 for large /x/ in M. Abramowitz and I.A. Stegun,
128 ['Handbook of Mathematical Functions] 9.2.19 are used
129 (these were found to be more reliable
130 than those in A&S 9.2.5).
132 When the order /v/ is an integer the method first relates the result
133 to J[sub 0], J[sub 1], Y[sub 0][space] and Y[sub 1][space] using either
134 forwards or backwards recurrence (Miller's algorithm) depending upon which is stable.
135 The values for J[sub 0], J[sub 1], Y[sub 0][space] and Y[sub 1][space] are
136 calculated using the rational minimax approximations on
137 root-bracketing intervals for small ['|x|] and Hankel asymptotic
138 expansion for large ['|x|]. The coefficients are from:
140 W.J. Cody, ['ALGORITHM 715: SPECFUN - A Portable FORTRAN Package of
141 Special Function Routines and Test Drivers], ACM Transactions on Mathematical
142 Software, vol 19, 22 (1993).
146 J.F. Hart et al, ['Computer Approximations], John Wiley & Sons, New York, 1968.
148 These approximations are accurate to around 19 decimal digits: therefore
149 these methods are not used when type T has more than 64 binary digits.
151 When /x/ is smaller than machine epsilon then the following approximations for
152 Y[sub 0](x), Y[sub 1](x), Y[sub 2](x) and Y[sub n](x) can be used
153 (see: [@http://functions.wolfram.com/03.03.06.0037.01 http://functions.wolfram.com/03.03.06.0037.01],
154 [@http://functions.wolfram.com/03.03.06.0038.01 http://functions.wolfram.com/03.03.06.0038.01],
155 [@http://functions.wolfram.com/03.03.06.0039.01 http://functions.wolfram.com/03.03.06.0039.01]
156 and [@http://functions.wolfram.com/03.03.06.0040.01 http://functions.wolfram.com/03.03.06.0040.01]):
158 [equation bessel_y0_small_z]
160 [equation bessel_y1_small_z]
162 [equation bessel_y2_small_z]
164 [equation bessel_yn_small_z]
166 When /x/ is small compared to /v/ and /v/ is not an integer, then the following
167 series approximation can be used for Y[sub v](x), this is also an area where other
168 approximations are often too slow to converge to be used
169 (see [@http://functions.wolfram.com/03.03.06.0034.01 http://functions.wolfram.com/03.03.06.0034.01]):
171 [equation bessel_yv_small_z]
173 When /x/ is small compared to /v/, J[sub v]x[space] is best computed directly from the series:
177 In the general case we compute J[sub v][space] and
178 Y[sub v][space] simultaneously.
180 To get the initial values, let
181 [mu][space] = [nu] - floor([nu] + 1/2), then [mu][space] is the fractional part
182 of [nu][space] such that
183 |[mu]| <= 1/2 (we need this for convergence later). The idea is to
184 calculate J[sub [mu]](x), J[sub [mu]+1](x), Y[sub [mu]](x), Y[sub [mu]+1](x)
185 and use them to obtain J[sub [nu]](x), Y[sub [nu]](x).
187 The algorithm is called Steed's method, which needs two
188 continued fractions as well as the Wronskian:
196 See: F.S. Acton, ['Numerical Methods that Work],
197 The Mathematical Association of America, Washington, 1997.
199 The continued fractions are computed using the modified Lentz's method
200 (W.J. Lentz, ['Generating Bessel functions in Mie scattering calculations
201 using continued fractions], Applied Optics, vol 15, 668 (1976)).
202 Their convergence rates depend on ['x], therefore we need
203 different strategies for large ['x] and small ['x].
205 ['x > v], CF1 needs O(['x]) iterations to converge, CF2 converges rapidly
207 ['x <= v], CF1 converges rapidly, CF2 fails to converge when ['x] [^->] 0
209 When ['x] is large (['x] > 2), both continued fractions converge (CF1
210 may be slow for really large ['x]). J[sub [mu]], J[sub [mu]+1],
211 Y[sub [mu]], Y[sub [mu]+1] can be calculated by
219 J[sub [nu]] and Y[sub [mu]] are then calculated using backward
220 (Miller's algorithm) and forward recurrence respectively.
222 When ['x] is small (['x] <= 2), CF2 convergence may fail (but CF1
223 works very well). The solution here is Temme's series:
231 g[sub k][space] and h[sub k][space]
232 are also computed by recursions (involving gamma functions), but the
233 formulas are a little complicated, readers are refered to
234 N.M. Temme, ['On the numerical evaluation of the ordinary Bessel function
235 of the second kind], Journal of Computational Physics, vol 21, 343 (1976).
236 Note Temme's series converge only for |[mu]| <= 1/2.
238 As the previous case, Y[sub [nu]][space] is calculated from the forward
239 recurrence, so is Y[sub [nu]+1]. With these two
240 values and f[sub [nu]], the Wronskian yields J[sub [nu]](x) directly
241 without backward recurrence.
245 [section:bessel_root Finding Zeros of Bessel Functions of the First and Second Kinds]
249 `#include <boost/math/special_functions/bessel.hpp>`
251 Functions for obtaining both a single zero or root of the Bessel function,
252 and placing multiple zeros into a container like `std::vector`
253 by providing an output iterator.
255 The signature of the single value functions are:
259 T v, // Floating-point value for Jv.
260 int m); // 1-based index of zero.
264 T v, // Floating-point value for Jv.
265 int m); // 1-based index of zero.
267 and for multiple zeros:
269 template <class T, class OutputIterator>
270 OutputIterator cyl_bessel_j_zero(
271 T v, // Floating-point value for Jv.
272 int start_index, // 1-based index of first zero.
273 unsigned number_of_zeros, // How many zeros to generate.
274 OutputIterator out_it); // Destination for zeros.
276 template <class T, class OutputIterator>
277 OutputIterator cyl_neumann_zero(
278 T v, // Floating-point value for Jv.
279 int start_index, // 1-based index of zero.
280 unsigned number_of_zeros, // How many zeros to generate
281 OutputIterator out_it); // Destination for zeros.
283 There are also versions which allow control of the __policy_section for error handling and precision.
287 T v, // Floating-point value for Jv.
288 int m, // 1-based index of zero.
289 const Policy&); // Policy to use.
293 T v, // Floating-point value for Jv.
294 int m, // 1-based index of zero.
295 const Policy&); // Policy to use.
298 template <class T, class OutputIterator>
299 OutputIterator cyl_bessel_j_zero(
300 T v, // Floating-point value for Jv.
301 int start_index, // 1-based index of first zero.
302 unsigned number_of_zeros, // How many zeros to generate.
303 OutputIterator out_it, // Destination for zeros.
304 const Policy& pol); // Policy to use.
306 template <class T, class OutputIterator>
307 OutputIterator cyl_neumann_zero(
308 T v, // Floating-point value for Jv.
309 int start_index, // 1-based index of zero.
310 unsigned number_of_zeros, // How many zeros to generate.
311 OutputIterator out_it, // Destination for zeros.
312 const Policy& pol); // Policy to use.
316 Every real order [nu] cylindrical Bessel and Neumann functions have an infinite
317 number of zeros on the positive real axis. The real zeros on the positive real
318 axis can be found by solving for the roots of
320 [emquad] ['J[sub [nu]](j[sub [nu], m]) = 0]
322 [emquad] ['Y[sub [nu]](y[sub [nu], m]) = 0]
324 Here, ['j[sub [nu], m]] represents the ['m[super th]]
325 root of the cylindrical Bessel function of order ['[nu]],
326 and ['y[sub [nu], m]] represents the ['m[super th]]
327 root of the cylindrical Neumann function of order ['[nu]].
329 The zeros or roots (values of `x` where the function crosses the horizontal `y = 0` axis)
330 of the Bessel and Neumann functions are computed by two functions,
331 `cyl_bessel_j_zero` and `cyl_neumann_zero`.
333 In each case the index or rank of the zero
334 returned is 1-based, which is to say:
336 cyl_bessel_j_zero(v, 1);
338 returns the first zero of Bessel J.
340 Passing an `start_index <= 0` results in a `std::domain_error` being raised.
342 For certain parameters, however, the zero'th root is defined and
343 it has a value of zero. For example, the zero'th root
344 of `J[sub v](x)` is defined and it has a value of zero for all
345 values of `v > 0` and for negative integer values of `v = -n`.
346 Similar cases are described in the implementation details below.
348 The order `v` of `J` can be positive, negative and zero for the `cyl_bessel_j`
349 and `cyl_neumann` functions, but not infinite nor NaN.
351 [graph bessel_j_zeros]
353 [graph neumann_y_zeros]
355 [h4 Examples of finding Bessel and Neumann zeros]
357 [import ../../example/bessel_zeros_example_1.cpp]
359 [bessel_zeros_example_1]
360 [bessel_zeros_example_2]
362 [import ../../example/bessel_zeros_interator_example.cpp]
364 [bessel_zeros_iterator_example_1]
365 [bessel_zeros_iterator_example_2]
367 [import ../../example/neumann_zeros_example_1.cpp]
369 [neumann_zeros_example_1]
370 [neumann_zeros_example_2]
372 [import ../../example/bessel_errors_example.cpp]
374 [bessel_errors_example_1]
375 [bessel_errors_example_2]
377 The full code (and output) for these examples is at
378 [@../../example/bessel_zeros_example_1.cpp Bessel zeros],
379 [@../../example/bessel_zeros_interator_example.cpp Bessel zeros iterator],
380 [@../../example/neumann_zeros_example_1.cpp Neumann zeros],
381 [@../../example/bessel_errors_example.cpp Bessel error messages].
385 Various methods are used to compute initial estimates
386 for ['j[sub [nu], m]] and ['y[sub [nu], m]] ; these are described in detail below.
388 After finding the initial estimate of a given root,
389 its precision is subsequently refined to the desired level
390 using Newton-Raphson iteration from Boost.Math's __root_finding_with_derivatives
391 utilities combined with the functions __cyl_bessel_j and __cyl_neumann.
393 Newton iteration requires both ['J[sub [nu]](x)] or ['Y[sub [nu]](x)]
394 as well as its derivative. The derivatives of ['J[sub [nu]](x)] and ['Y[sub [nu]](x)]
395 with respect to ['x] are given by __Abramowitz_Stegun.
398 [emquad] ['d/[sub dx] ['J[sub [nu]](x)] = ['J[sub [nu]-1](x)] - [nu] J[sub [nu]](x)] / x
400 [emquad] ['d/[sub dx] ['Y[sub [nu]](x)] = ['Y[sub [nu]-1](x)] - [nu] Y[sub [nu]](x)] / x
402 Enumeration of the rank of a root (in other words the index of a root)
403 begins with one and counts up, in other words
404 ['m,=1,2,3,[ellipsis]] The value of the first root is always greater than zero.
406 For certain special parameters, cylindrical Bessel functions
407 and cylindrical Neumann functions have a root at the origin. For example,
408 ['J[sub [nu]](x)] has a root at the origin for every positive order
409 ['[nu] > 0], and for every negative integer order
410 ['[nu] = -n] with ['n [isin] [negative] [super +]] and ['n [ne] 0].
412 In addition, ['Y[sub [nu]](x)] has a root at the origin
413 for every negative half-integer order ['[nu] = -n/2], with ['n [isin] [negative] [super +]] and
416 For these special parameter values, the origin with
417 a value of ['x = 0] is provided as the ['0[super th]]
418 root generated by `cyl_bessel_j_zero()`
419 and `cyl_neumann_zero()`.
421 When calculating initial estimates for the roots
422 of Bessel functions, a distinction is made between
423 positive order and negative order, and different
424 methods are used for these. In addition, different algorithms
425 are used for the first root ['m = 1] and
426 for subsequent roots with higher rank ['m [ge] 2].
427 Furthermore, estimates of the roots for Bessel functions
428 with order above and below a cutoff at ['[nu] = 2.2]
429 are calculated with different methods.
431 Calculations of the estimates of ['j[sub [nu],1]] and ['y[sub [nu],1]]
432 with ['0 [le] [nu] < 2.2] use empirically tabulated values.
433 The coefficients for these have been generated by a
434 computer algebra system.
436 Calculations of the estimates of ['j[sub [nu],1]] and ['y[sub [nu],1]]
437 with ['[nu][ge] 2.2] use Eqs.9.5.14 and 9.5.15 in __Abramowitz_Stegun.
441 [emquad] ['j[sub [nu],1] [cong] [nu] + 1.85575 [nu][super [frac13]] + 1.033150 [nu][super -[frac13]] - 0.00397 [nu][super -1] - 0.0908 [nu][super -5/3] + 0.043 [nu][super -7/3] + [ellipsis]]
445 [emquad] ['y[sub [nu],1] [cong] [nu] + 0.93157 [nu][super [frac13]] + 0.26035 [nu][super -[frac13]] + 0.01198 [nu][super -1] - 0.0060 [nu][super -5/3] - 0.001 [nu][super -7/3] + [ellipsis]]
447 Calculations of the estimates of ['j[sub [nu], m]] and ['y[sub [nu], m]]
448 with rank ['m > 2] and ['0 [le] [nu] < 2.2] use
449 McMahon's approximation, as described in M. Abramowitz and I. A. Stegan, Section 9.5 and 9.5.12.
452 [emquad] ['j[sub [nu],m], y[sub [nu],m] [cong] [beta] - ([mu]-1) / 8[beta]]
454 [emquad] [emquad] [emquad] ['- 4([mu]-1)(7[mu] - 31) / 3(8[beta])[super 3]]
456 [emquad] [emquad] [emquad] ['-32([mu]-1)(83[mu][sup2] - 982[mu] + 3779) / 15(8[beta])[super 5]]
458 [emquad] [emquad] [emquad] ['-64([mu]-1)(6949[mu][super 3] - 153855[mu][sup2] + 1585743[mu]- 6277237) / 105(8a)[super 7]]
460 [emquad] [emquad] [emquad] ['- [ellipsis]] [emquad] [emquad] (5)
462 where ['[mu] = 4[nu][super 2]] and ['[beta] = (m + [frac12][nu] - [frac14])[pi]]
463 for ['j[sub [nu],m]] and
464 ['[beta] = (m + [frac12][nu] -[frac34])[pi] for ['y[sub [nu],m]]].
466 Calculations of the estimates of ['j[sub [nu], m]] and ['y[sub [nu], m]]
467 with ['[nu] [ge] 2.2] use
468 one term in the asymptotic expansion given in
469 Eq.9.5.22 and top line of Eq.9.5.26 combined with Eq. 9.3.39,
470 all in __Abramowitz_Stegun explicit and easy-to-understand treatment
471 for asymptotic expansion of zeros.
472 The latter two equations are expressed for argument ['(x)] greater than one.
473 (Olver also gives the series form of the equations in
474 [@http://dlmf.nist.gov/10.21#vi [sect]10.21(vi) McMahon's Asymptotic Expansions for Large Zeros] - using slightly different variable names).
478 [emquad] ['j[sub [nu], m] [sim] [nu]x(-[zeta]) + f[sub 1](-[zeta]/[nu])]
480 where ['-[zeta] = [nu][super -2/3]a[sub m]] and ['a[sub m]] is
481 the absolute value of the ['m[super th]] root of ['Ai(x)] on the negative real axis.
483 Here ['x = x(-[zeta])] is the inverse of the function
485 [emquad] ['[frac23](-[zeta])[super 3/2] = [radic](x[sup2] - 1) - cos[supminus][sup1](1/x)] [emquad] [emquad] (7)
489 [emquad] ['f[sub 1](-[zeta]) = [frac12]x(-[zeta]) {h(-[zeta])}[sup2] [sdot] b[sub 0](-[zeta])]
493 [emquad] ['h(-[zeta]) = {4(-[zeta]) / (x[sup2] - 1)}[super 4]]
497 [emquad] ['b[sub 0](-[zeta]) = -5/(48[zeta][sup2]) + 1/(-[zeta])[super [frac12]] [sdot] { 5/(24(x[super 2]-1)[super 3/2]) + 1/(8(x[super 2]-1)[super [frac12])]}]
499 When solving for ['x(-[zeta])] in Eq. 7 above,
500 the right-hand-side is expanded to order 2 in
501 a Taylor series for large ['x]. This results in
503 [emquad] ['[frac23](-[zeta])[super 3/2] [approx] x + 1/2x - [pi]/2]
505 The positive root of the resulting quadratic equation
506 is used to find an initial estimate ['x(-[zeta])].
507 This initial estimate is subsequently refined with
508 several steps of Newton-Raphson iteration
511 Estimates of the roots of cylindrical Bessel functions
512 of negative order on the positive real axis are found
513 using interlacing relations. For example, the ['m[super th]]
514 root of the cylindrical Bessel function ['j[sub -[nu],m]]
515 is bracketed by the ['m[super th]] root and the
516 ['(m+1)[super th]] root of the Bessel function of
517 corresponding positive integer order. In other words,
519 [emquad] ['j[sub n[nu],m]] < ['j[sub -[nu],m]] < ['j[sub n[nu],m+1]]
521 where ['m > 1] and ['n[sub [nu]]] represents the integral
522 floor of the absolute value of ['|-[nu]|].
524 Similar bracketing relations are used to find estimates
525 of the roots of Neumann functions of negative order,
526 whereby a discontinuity at every negative half-integer
527 order needs to be handled.
529 Bracketing relations do not hold for the first root
530 of cylindrical Bessel functions and cylindrical Neumann
531 functions with negative order. Therefore, iterative algorithms
532 combined with root-finding via bisection are used
533 to localize ['j[sub -[nu],1]] and ['y[sub -[nu],1]].
537 The precision of evaluation of zeros was tested at 50 decimal digits using `cpp_dec_float_50`
538 and found identical with spot values computed by __WolframAlpha.
540 [endsect] [/section:bessel Finding Zeros of Bessel Functions of the First and Second Kinds]
543 Copyright 2006, 2013 John Maddock, Paul A. Bristow, Xiaogang Zhang and Christopher Kormanyos.
545 Distributed under the Boost Software License, Version 1.0.
546 (See accompanying file LICENSE_1_0.txt or copy at
547 http://www.boost.org/LICENSE_1_0.txt).