libs/math/doc/sf/bessel_jy.qbk

   1 [section:bessel_first Bessel Functions of the First and Second Kinds]
   2
   3 [h4 Synopsis]
   4
   5 `#include <boost/math/special_functions/bessel.hpp>`
   6
   7    template <class T1, class T2>
   8    ``__sf_result`` cyl_bessel_j(T1 v, T2 x);
   9
  10    template <class T1, class T2, class ``__Policy``>
  11    ``__sf_result`` cyl_bessel_j(T1 v, T2 x, const ``__Policy``&);
  12
  13    template <class T1, class T2>
  14    ``__sf_result`` cyl_neumann(T1 v, T2 x);
  15
  16    template <class T1, class T2, class ``__Policy``>
  17    ``__sf_result`` cyl_neumann(T1 v, T2 x, const ``__Policy``&);
  18
  19
  20 [h4 Description]
  21
  22 The functions __cyl_bessel_j and __cyl_neumann return the result of the
  23 Bessel functions of the first and second kinds respectively:
  24
  25 [expression cyl_bessel_j(v, x) = J[sub v](x)]
  26
  27 [expression cyl_neumann(v, x) = Y[sub v](x) = N[sub v](x)]
  28
  29 where:
  30
  31 [equation bessel2]
  32
  33 [equation bessel3]
  34
  35 The return type of these functions is computed using the __arg_promotion_rules
  36 when T1 and T2 are different types.  The functions are also optimised for the
  37 relatively common case that T1 is an integer.
  38
  39 [optional_policy]
  40
  41 The functions return the result of __domain_error whenever the result is
  42 undefined or complex.  For __cyl_bessel_j this occurs when `x < 0` and v is not
  43 an integer, or when `x == 0` and `v != 0`.  For __cyl_neumann this occurs
  44 when `x <= 0`.
  45
  46 The following graph illustrates the cyclic nature of J[sub v]:
  47
  48 [graph cyl_bessel_j]
  49
  50 The following graph shows the behaviour of Y[sub v]: this is also
  51 cyclic for large /x/, but tends to -[infin] for small /x/:
  52
  53 [graph cyl_neumann]
  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 comparisons to other
  67 libraries.  Note that the cyclic nature of these
  68 functions means that they have an infinite number of irrational
  69 roots: in general these functions have arbitrarily large /relative/
  70 errors when the arguments are sufficiently close to a root.  Of
  71 course the absolute error in such cases is always small.
  72 Note that only results for the widest floating-point type on the
  73 system are given as narrower types have __zero_error.  All values
  74 are relative errors in units of epsilon.  Most of the gross errors
  75 exhibited by other libraries occur for very large arguments - you will
  76 need to drill down into the actual program output if you need more
  77 information on this.
  78
  79 [table_cyl_bessel_j_integer_orders_]
  80
  81 [table_cyl_bessel_j]
  82
  83 [table_cyl_neumann_integer_orders_]
  84
  85 [table_cyl_neumann]
  86
  87 Note that for large /x/ these functions are largely dependent on
  88 the accuracy of the `std::sin` and `std::cos` functions.
  89
  90 Comparison to GSL and __cephes is interesting: both __cephes and this library optimise
  91 the integer order case - leading to identical results - simply using the general
  92 case is for the most part slightly more accurate though, as noted by the
  93 better accuracy of GSL in the integer argument cases.  This implementation tends to
  94 perform much better when the arguments become large, __cephes in particular produces
  95 some remarkably inaccurate results with some of the test data (no significant figures
  96 correct), and even GSL performs badly with some inputs to J[sub v].  Note that
  97 by way of double-checking these results, the worst performing __cephes and GSL cases
  98 were recomputed using [@http://functions.wolfram.com functions.wolfram.com],
  99 and the result checked against our test data: no errors in the test data were found.
 100
 101 The following error plot are based on an exhaustive search of the functions domain for J0 and Y0,
 102 MSVC-15.5 at `double` precision, other compilers and precisions are very similar - the plots simply
 103 illustrate the relatively large errors as you approach a zero, and the very low errors elsewhere.
 104
 105 [graph j0__double]
 106
 107 [graph y0__double]
 108
 109
 110 [h4 Implementation]
 111
 112 The implementation is mostly about filtering off various special cases:
 113
 114 When /x/ is negative, then the order /v/ must be an integer or the
 115 result is a domain error.  If the order is an integer then the function
 116 is odd for odd orders and even for even orders, so we reflect to /x > 0/.
 117
 118 When the order /v/ is negative then the reflection formulae can be used to
 119 move to /v > 0/:
 120
 121 [equation bessel9]
 122
 123 [equation bessel10]
 124
 125 Note that if the order is an integer, then these formulae reduce to:
 126
 127 [expression J[sub -n] = (-1)[super n]J[sub n]]
 128
 129 [expression Y[sub -n] = (-1)[super n]Y[sub n]]
 130
 131 However, in general, a negative order implies that we will need to compute
 132 both J and Y.
 133
 134 When /x/ is large compared to the order /v/ then the asymptotic expansions
 135 for large /x/ in M. Abramowitz and I.A. Stegun,
 136 ['Handbook of Mathematical Functions] 9.2.19 are used
 137 (these were found to be more reliable
 138 than those in A&S 9.2.5).
 139
 140 When the order /v/ is an integer the method first relates the result
 141 to J[sub 0], J[sub 1], Y[sub 0] and Y[sub 1] using either
 142 forwards or backwards recurrence (Miller's algorithm) depending upon which is stable.
 143 The values for J[sub 0], J[sub 1], Y[sub 0] and Y[sub 1] are
 144 calculated using the rational minimax approximations on
 145 root-bracketing intervals for small ['|x|] and Hankel asymptotic
 146 expansion for large ['|x|]. The coefficients are from:
 147
 148 [:W.J. Cody, ['ALGORITHM 715: SPECFUN - A Portable FORTRAN Package of
 149 Special Function Routines and Test Drivers], ACM Transactions on Mathematical
 150 Software, vol 19, 22 (1993).]
 151
 152 and
 153
 154 [:J.F. Hart et al, ['Computer Approximations], John Wiley & Sons, New York, 1968.]
 155
 156 These approximations are accurate to around 19 decimal digits: therefore
 157 these methods are not used when type T has more than 64 binary digits.
 158
 159 When /x/ is smaller than machine epsilon then the following approximations for
 160 Y[sub 0](x), Y[sub 1](x), Y[sub 2](x) and Y[sub n](x) can be used
 161 (see: [@http://functions.wolfram.com/03.03.06.0037.01 http://functions.wolfram.com/03.03.06.0037.01],
 162 [@http://functions.wolfram.com/03.03.06.0038.01 http://functions.wolfram.com/03.03.06.0038.01],
 163 [@http://functions.wolfram.com/03.03.06.0039.01 http://functions.wolfram.com/03.03.06.0039.01]
 164 and [@http://functions.wolfram.com/03.03.06.0040.01 http://functions.wolfram.com/03.03.06.0040.01]):
 165
 166 [equation bessel_y0_small_z]
 167
 168 [equation bessel_y1_small_z]
 169
 170 [equation bessel_y2_small_z]
 171
 172 [equation bessel_yn_small_z]
 173
 174 When /x/ is small compared to /v/ and /v/ is not an integer, then the following
 175 series approximation can be used for Y[sub v](x), this is also an area where other
 176 approximations are often too slow to converge to be used
 177 (see [@http://functions.wolfram.com/03.03.06.0034.01 http://functions.wolfram.com/03.03.06.0034.01]):
 178
 179 [equation bessel_yv_small_z]
 180
 181 When /x/ is small compared to /v/, J[sub v]x is best computed directly from the series:
 182
 183 [equation bessel2]
 184
 185 In the general case we compute J[sub v] and
 186 Y[sub v] simultaneously.
 187
 188 To get the initial values, let
 189 [mu] = [nu] - floor([nu] + 1/2), then [mu] is the fractional part
 190 of [nu] such that
 191 |[mu]| <= 1/2 (we need this for convergence later). The idea is to
 192 calculate J[sub [mu]](x), J[sub [mu]+1](x), Y[sub [mu]](x), Y[sub [mu]+1](x)
 193 and use them to obtain J[sub [nu]](x), Y[sub [nu]](x).
 194
 195 The algorithm is called Steed's method, which needs two
 196 continued fractions as well as the Wronskian:
 197
 198 [equation bessel8]
 199
 200 [equation bessel11]
 201
 202 [equation bessel12]
 203
 204 See: F.S. Acton, ['Numerical Methods that Work],
 205     The Mathematical Association of America, Washington, 1997.
 206
 207 The continued fractions are computed using the modified Lentz's method
 208 (W.J. Lentz, ['Generating Bessel functions in Mie scattering calculations
 209 using continued fractions], Applied Optics, vol 15, 668 (1976)).
 210 Their convergence rates depend on ['x], therefore we need
 211 different strategies for large ['x] and small ['x]:
 212
 213 [:['x > v], CF1 needs O(['x]) iterations to converge, CF2 converges rapidly]
 214
 215 [:['x <= v], CF1 converges rapidly, CF2 fails to converge when ['x] [^->] 0]
 216
 217 When ['x] is large (['x] > 2), both continued fractions converge (CF1
 218 may be slow for really large ['x]). J[sub [mu]], J[sub [mu]+1],
 219 Y[sub [mu]], Y[sub [mu]+1] can be calculated by
 220
 221 [equation bessel13]
 222
 223 where
 224
 225 [equation bessel14]
 226
 227 J[sub [nu]] and Y[sub [mu]] are then calculated using backward
 228 (Miller's algorithm) and forward recurrence respectively.
 229
 230 When ['x] is small (['x] <= 2), CF2 convergence may fail (but CF1
 231 works very well). The solution here is Temme's series:
 232
 233 [equation bessel15]
 234
 235 where
 236
 237 [equation bessel16]
 238
 239 g[sub k] and h[sub k]
 240 are also computed by recursions (involving gamma functions), but the
 241 formulas are a little complicated, readers are refered to
 242 N.M. Temme, ['On the numerical evaluation of the ordinary Bessel function
 243 of the second kind], Journal of Computational Physics, vol 21, 343 (1976).
 244 Note Temme's series converge only for |[mu]| <= 1/2.
 245
 246 As the previous case, Y[sub [nu]] is calculated from the forward
 247 recurrence, so is Y[sub [nu]+1]. With these two
 248 values and f[sub [nu]], the Wronskian yields J[sub [nu]](x) directly
 249 without backward recurrence.
 250
 251 [endsect] [/section:bessel_first Bessel Functions of the First and Second Kinds]
 252
 253 [section:bessel_root Finding Zeros of Bessel Functions of the First and Second Kinds]
 254
 255 [h4 Synopsis]
 256
 257 `#include <boost/math/special_functions/bessel.hpp>`
 258
 259 Functions for obtaining both a single zero or root of the Bessel function,
 260 and placing multiple zeros into a container like `std::vector`
 261 by providing an output iterator.
 262
 263 The signature of the single value functions are:
 264
 265   template <class T>
 266   T cyl_bessel_j_zero(
 267            T v,            // Floating-point value for Jv.
 268            int m);         // 1-based index of zero.
 269
 270   template <class T>
 271   T cyl_neumann_zero(
 272            T v,            // Floating-point value for Jv.
 273            int m);         // 1-based index of zero.
 274
 275 and for multiple zeros:
 276
 277  template <class T, class OutputIterator>
 278  OutputIterator cyl_bessel_j_zero(
 279                       T v,                       // Floating-point value for Jv.
 280                       int start_index,           // 1-based index of first zero.
 281                       unsigned number_of_zeros,  // How many zeros to generate.
 282                       OutputIterator out_it);    // Destination for zeros.
 283
 284  template <class T, class OutputIterator>
 285  OutputIterator cyl_neumann_zero(
 286                       T v,                       // Floating-point value for Jv.
 287                       int start_index,           // 1-based index of zero.
 288                       unsigned number_of_zeros,  // How many zeros to generate
 289                       OutputIterator out_it);    // Destination for zeros.
 290
 291 There are also versions which allow control of the __policy_section for error handling and precision.
 292
 293   template <class T>
 294   T cyl_bessel_j_zero(
 295            T v,            // Floating-point value for Jv.
 296            int m,          // 1-based index of zero.
 297            const Policy&); // Policy to use.
 298
 299   template <class T>
 300   T cyl_neumann_zero(
 301            T v,            // Floating-point value for Jv.
 302            int m,          // 1-based index of zero.
 303            const Policy&); // Policy to use.
 304
 305
 306  template <class T, class OutputIterator>
 307  OutputIterator cyl_bessel_j_zero(
 308                       T v,                       // Floating-point value for Jv.
 309                       int start_index,           // 1-based index of first 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.
 313
 314  template <class T, class OutputIterator>
 315  OutputIterator cyl_neumann_zero(
 316                       T v,                       // Floating-point value for Jv.
 317                       int start_index,           // 1-based index of zero.
 318                       unsigned number_of_zeros,  // How many zeros to generate.
 319                       OutputIterator out_it,     // Destination for zeros.
 320                       const Policy& pol);        // Policy to use.
 321
 322 [h4 Description]
 323
 324 Every real order [nu] cylindrical Bessel and Neumann functions have an infinite
 325 number of zeros on the positive real axis. The real zeros on the positive real
 326 axis can be found by solving for the roots of
 327
 328 [:['J[sub [nu]](j[sub [nu], m]) = 0]]
 329
 330 [:['Y[sub [nu]](y[sub [nu], m]) = 0]]
 331
 332 Here, ['j[sub [nu], m]] represents the ['m[super th]]
 333 root of the cylindrical Bessel function of order ['[nu]],
 334 and ['y[sub [nu], m]] represents the ['m[super th]]
 335 root of the cylindrical Neumann function of order ['[nu]].
 336
 337 The zeros or roots (values of `x` where the function crosses the horizontal `y = 0` axis)
 338 of the Bessel and Neumann functions are computed by two functions,
 339 `cyl_bessel_j_zero` and `cyl_neumann_zero`.
 340
 341 In each case the index or rank of the zero
 342 returned is 1-based, which is to say:
 343
 344 [:cyl_bessel_j_zero(v, 1);]
 345
 346 returns the first zero of Bessel J.
 347
 348 Passing an `start_index <= 0` results in a `std::domain_error` being raised.
 349
 350 For certain parameters, however, the zero'th root is defined and
 351 it has a value of zero. For example, the zero'th root
 352 of `J[sub v](x)` is defined and it has a value of zero for all
 353 values of `v > 0` and for negative integer values of `v = -n`.
 354 Similar cases are described in the implementation details below.
 355
 356 The order `v` of `J` can be positive, negative and zero for the `cyl_bessel_j`
 357 and `cyl_neumann` functions, but not infinite nor NaN.
 358
 359 [graph bessel_j_zeros]
 360
 361 [graph neumann_y_zeros]
 362
 363 [h4 Examples of finding Bessel and Neumann zeros]
 364
 365 [import ../../example/bessel_zeros_example_1.cpp]
 366
 367 [bessel_zeros_example_1]
 368 [bessel_zeros_example_2]
 369
 370 [import ../../example/bessel_zeros_interator_example.cpp]
 371
 372 [bessel_zeros_iterator_example_1]
 373 [bessel_zeros_iterator_example_2]
 374
 375 [import ../../example/neumann_zeros_example_1.cpp]
 376
 377 [neumann_zeros_example_1]
 378 [neumann_zeros_example_2]
 379
 380 [import ../../example/bessel_errors_example.cpp]
 381
 382 [bessel_errors_example_1]
 383 [bessel_errors_example_2]
 384
 385 The full code (and output) for these examples is at
 386 [@../../example/bessel_zeros_example_1.cpp Bessel zeros],
 387 [@../../example/bessel_zeros_interator_example.cpp Bessel zeros iterator],
 388 [@../../example/neumann_zeros_example_1.cpp Neumann zeros],
 389 [@../../example/bessel_errors_example.cpp  Bessel error messages].
 390
 391 [h3 Implementation]
 392
 393 Various methods are used to compute initial estimates
 394 for ['j[sub [nu], m]] and ['y[sub [nu], m]] ; these are described in detail below.
 395
 396 After finding the initial estimate of a given root,
 397 its precision is subsequently refined to the desired level
 398 using Newton-Raphson iteration from Boost.Math's __root_finding_with_derivatives
 399 utilities combined with the functions __cyl_bessel_j  and __cyl_neumann.
 400
 401 Newton iteration requires both ['J[sub [nu]](x)] or ['Y[sub [nu]](x)]
 402 as well as its derivative. The derivatives of ['J[sub [nu]](x)] and ['Y[sub [nu]](x)]
 403 with respect to  ['x] are given by __Abramowitz_Stegun.
 404 In particular,
 405
 406 [expression d/[sub dx] ['J[sub [nu]](x)] = ['J[sub [nu]-1](x)] - [nu] J[sub [nu]](x) / x]
 407
 408 [expression d/[sub dx] ['Y[sub [nu]](x)] = ['Y[sub [nu]-1](x)] - [nu] Y[sub [nu]](x) / x]
 409
 410 Enumeration of the rank of a root (in other words the index of a root)
 411 begins with one and counts up, in other words
 412 ['m,=1,2,3,[ellipsis]] The value of the first root is always greater than zero.
 413
 414 For certain special parameters, cylindrical Bessel functions
 415 and cylindrical Neumann functions have a root at the origin. For example,
 416 ['J[sub [nu]](x)] has a root at the origin for every positive order
 417 ['[nu] > 0], and for every negative integer order
 418 ['[nu] = -n] with ['n [isin] [negative] [super +]] and ['n [ne] 0].
 419
 420 In addition, ['Y[sub [nu]](x)] has a root at the origin
 421 for every negative half-integer order ['[nu] = -n/2], with ['n [isin] [negative] [super +]] and
 422 and ['n [ne] 0].
 423
 424 For these special parameter values, the origin with
 425 a value of ['x = 0] is provided as the ['0[super th]]
 426 root generated by `cyl_bessel_j_zero()`
 427 and `cyl_neumann_zero()`.
 428
 429 When calculating initial estimates for the roots
 430 of Bessel functions, a distinction is made between
 431 positive order and negative order, and different
 432 methods are used for these. In addition, different algorithms
 433 are used for the first root ['m = 1] and
 434 for subsequent roots with higher rank ['m [ge] 2].
 435 Furthermore, estimates of the roots for Bessel functions
 436 with order above and below a cutoff at ['[nu] = 2.2]
 437 are calculated with different methods.
 438
 439 Calculations of the estimates of ['j[sub [nu],1]] and ['y[sub [nu],1]]
 440 with ['0 [le] [nu] < 2.2] use empirically tabulated values.
 441 The coefficients for these have been generated by a
 442 computer algebra system.
 443
 444 Calculations of the estimates of ['j[sub [nu],1]] and ['y[sub [nu],1]]
 445 with ['[nu][ge] 2.2] use Eqs.9.5.14 and 9.5.15 in __Abramowitz_Stegun.
 446
 447 In particular,
 448 [expression 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]]
 449
 450 and
 451
 452 [expression 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]]
 453
 454 Calculations of the estimates of ['j[sub [nu], m]]  and  ['y[sub [nu], m]]
 455 with rank ['m > 2] and ['0 [le] [nu] < 2.2]  use
 456 McMahon's approximation, as described in M. Abramowitz and I. A. Stegan, Section 9.5 and 9.5.12.
 457 In particular,
 458
 459 [:['j[sub [nu],m], y[sub [nu],m] [cong]]]
 460 [:[:[beta] - ([mu]-1) / 8[beta]]]
 461 [:[:['- 4([mu]-1)(7[mu] - 31) / 3(8[beta])[super 3]]]]
 462 [:[:['-32([mu]-1)(83[mu][sup2] - 982[mu] + 3779) / 15(8[beta])[super 5]]]]
 463 [:[:['-64([mu]-1)(6949[mu][super 3] - 153855[mu][sup2] + 1585743[mu]- 6277237) / 105(8a)[super 7]]]]
 464 [:[:['- [ellipsis]] [emquad] (5)]]
 465
 466 where ['[mu] = 4[nu][super 2]] and ['[beta] = (m + [frac12][nu] - [frac14])[pi]]
 467 for ['j[sub [nu],m]] and
 468 ['[beta] = (m + [frac12][nu] -[frac34])[pi] for ['y[sub [nu],m]]].
 469
 470 Calculations of the estimates of ['j[sub [nu], m]]  and  ['y[sub [nu], m]]
 471 with ['[nu] [ge] 2.2] use
 472 one term in the asymptotic expansion given in
 473 Eq.9.5.22 and top line of Eq.9.5.26 combined with Eq. 9.3.39,
 474 all in __Abramowitz_Stegun explicit and easy-to-understand treatment
 475 for asymptotic expansion of zeros.
 476 The latter two equations are expressed for argument ['(x)] greater than one.
 477 (Olver also gives the series form of the equations in
 478 [@http://dlmf.nist.gov/10.21#vi [sect]10.21(vi) McMahon's Asymptotic Expansions for Large Zeros] - using slightly different variable names).
 479
 480 In summary,
 481
 482 [expression j[sub [nu], m] [sim] [nu]x(-[zeta]) + f[sub 1](-[zeta]/[nu])]
 483
 484 where ['-[zeta] = [nu][super -2/3]a[sub m]] and ['a[sub m]] is
 485 the absolute value of the ['m[super th]] root of ['Ai(x)] on the negative real axis.
 486
 487 Here ['x = x(-[zeta])] is the inverse of the function
 488
 489 [expression [frac23](-[zeta])[super 3/2] = [radic](x[sup2] - 1) - cos[supminus][sup1](1/x)]     (7)
 490
 491 Furthermore,
 492
 493 [expression f[sub 1](-[zeta]) = [frac12]x(-[zeta]) {h(-[zeta])}[sup2] [sdot] b[sub 0](-[zeta])]
 494
 495 where
 496
 497 [expression h(-[zeta]) = {4(-[zeta]) / (x[sup2] - 1)}[super 4]]
 498
 499 and
 500
 501 [expression 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])]}]
 502
 503 When solving for ['x(-[zeta])] in Eq. 7 above,
 504 the right-hand-side is expanded to order 2 in
 505 a Taylor series for large ['x]. This results in
 506
 507 [expression [frac23](-[zeta])[super 3/2] [approx] x + 1/2x - [pi]/2]
 508
 509 The positive root of the resulting quadratic equation
 510 is used to find an initial estimate ['x(-[zeta])].
 511 This initial estimate is subsequently refined with
 512 several steps of Newton-Raphson iteration
 513 in Eq. 7.
 514
 515 Estimates of the roots of cylindrical Bessel functions
 516 of negative order on the positive real axis are found
 517 using interlacing relations. For example, the ['m[super th]]
 518 root of the cylindrical Bessel function ['j[sub -[nu],m]]
 519 is bracketed by the ['m[super th]] root and the
 520 ['(m+1)[super th]] root of the Bessel function of
 521 corresponding positive integer order. In other words,
 522
 523 [expression j[sub n[nu], m] < j[sub -[nu], m] < j[sub n[nu], m+1]]
 524
 525 where ['m > 1] and ['n[sub [nu]]] represents the integral
 526 floor of the absolute value of ['|-[nu]|].
 527
 528 Similar bracketing relations are used to find estimates
 529 of the roots of Neumann functions of negative order,
 530 whereby a discontinuity at every negative half-integer
 531 order needs to be handled.
 532
 533 Bracketing relations do not hold for the first root
 534 of cylindrical Bessel functions and cylindrical Neumann
 535 functions with negative order. Therefore, iterative algorithms
 536 combined with root-finding via bisection are used
 537 to localize ['j[sub -[nu],1]] and ['y[sub -[nu],1]].
 538
 539 [h3 Testing]
 540
 541 The precision of evaluation of zeros was tested at 50 decimal digits using `cpp_dec_float_50`
 542 and found identical with spot values computed by __WolframAlpha.
 543
 544 [endsect]  [/section:bessel Finding Zeros of Bessel Functions of the First and Second Kinds]
 545
 546 [/
 547   Copyright 2006, 2013 John Maddock, Paul A. Bristow, Xiaogang Zhang and Christopher Kormanyos.
 548
 549   Distributed under the Boost Software License, Version 1.0.
 550   (See accompanying file LICENSE_1_0.txt or copy at
 551   http://www.boost.org/LICENSE_1_0.txt).
 552 ]