| |
| [section:bessel Bessel Functions of the First and Second Kinds] |
| |
| [h4 Synopsis] |
| |
| template <class T1, class T2> |
| ``__sf_result`` cyl_bessel_j(T1 v, T2 x); |
| |
| template <class T1, class T2, class ``__Policy``> |
| ``__sf_result`` cyl_bessel_j(T1 v, T2 x, const ``__Policy``&); |
| |
| template <class T1, class T2> |
| ``__sf_result`` cyl_neumann(T1 v, T2 x); |
| |
| template <class T1, class T2, class ``__Policy``> |
| ``__sf_result`` cyl_neumann(T1 v, T2 x, const ``__Policy``&); |
| |
| |
| [h4 Description] |
| |
| The functions __cyl_bessel_j and __cyl_neumann return the result of the |
| Bessel functions of the first and second kinds respectively: |
| |
| cyl_bessel_j(v, x) = J[sub v](x) |
| |
| cyl_neumann(v, x) = Y[sub v](x) = N[sub v](x) |
| |
| where: |
| |
| [equation bessel2] |
| |
| [equation bessel3] |
| |
| The return type of these functions is computed using the __arg_pomotion_rules |
| when T1 and T2 are different types. The functions are also optimised for the |
| relatively common case that T1 is an integer. |
| |
| [optional_policy] |
| |
| The functions return the result of __domain_error whenever the result is |
| undefined or complex. For __cyl_bessel_j this occurs when `x < 0` and v is not |
| an integer, or when `x == 0` and `v != 0`. For __cyl_neumann this occurs |
| when `x <= 0`. |
| |
| The following graph illustrates the cyclic nature of J[sub v]: |
| |
| [graph cyl_bessel_j] |
| |
| The following graph shows the behaviour of Y[sub v]: this is also |
| cyclic for large /x/, but tends to -[infin][space] for small /x/: |
| |
| [graph cyl_neumann] |
| |
| [h4 Testing] |
| |
| There are two sets of test values: spot values calculated using |
| [@http://functions.wolfram.com functions.wolfram.com], |
| and a much larger set of tests computed using |
| a simplified version of this implementation |
| (with all the special case handling removed). |
| |
| [h4 Accuracy] |
| |
| The following tables show how the accuracy of these functions |
| varies on various platforms, along with comparisons to the __gsl and |
| __cephes 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 /relative/ |
| 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 __zero_error. All values |
| are relative errors in units of epsilon. |
| |
| [table Errors Rates in cyl_bessel_j |
| [[Significand Size] [Platform and Compiler] [J[sub 0][space] and J[sub 1]] [J[sub v]] [J[sub v][space] (large values of x > 1000)] ] |
| [[53] [Win32 / Visual C++ 8.0] |
| [Peak=2.5 Mean=1.1 |
| |
| GSL Peak=6.6 |
| |
| __cephes Peak=2.5 Mean=1.1] |
| [Peak=11 Mean=2.2 |
| |
| GSL Peak=11 |
| |
| __cephes Peak=17 Mean=2.5] |
| [Peak=413 Mean=110 |
| |
| GSL Peak=6x10[super 11] |
| |
| __cephes Peak=2x10[super 5] ] ] |
| [[64] [Red Hat Linux IA64 / G++ 3.4] [Peak=7 Mean=3] [Peak=117 Mean=10] [Peak=2x10[super 4][space] Mean=6x10[super 3]] ] |
| [[64] [SUSE Linux AMD64 / G++ 4.1] [Peak=7 Mean=3] [Peak=400 Mean=40] [Peak=2x10[super 4][space] Mean=1x10[super 4]] ] |
| [[113] [HP-UX / HP aCC 6] [Peak=14 Mean=6] [Peak=29 Mean=3] [Peak=2700 Mean=450] ] |
| ] |
| |
| [table Errors Rates in cyl_neumann |
| [[Significand Size] [Platform and Compiler] [J[sub 0][space] and J[sub 1]] [J[sub n] (integer orders)] [J[sub v] (fractional orders)] ] |
| [[53] [Win32 / Visual C++ 8.0] |
| [Peak=330 Mean=54 |
| |
| GSL Peak=34 Mean=9 |
| |
| __cephes Peak=330 Mean=54] |
| [Peak=923 Mean=83 |
| |
| GSL Peak=500 Mean=54 |
| |
| __cephes Peak=923 Mean=83] |
| [Peak=561 Mean=36 |
| |
| GSL Peak=1.4x10[super 6][space] Mean\=7x10[super 4][space] |
| |
| __cephes Peak=+INF]] |
| [[64] [Red Hat Linux IA64 / G++ 3.4] [Peak=470 Mean=56] [Peak=843 Mean=51] [Peak=741 Mean=51] ] |
| [[64] [SUSE Linux AMD64 / G++ 4.1] [Peak=1300 Mean=424] [Peak=2x10[super 4][space] Mean=8x10[super 3]] [Peak=1x10[super 5][space] Mean=6x10[super 3]] ] |
| [[113] [HP-UX / HP aCC 6] [Peak=180 Mean=63] [Peak=340 Mean=150] [Peak=2x10[super 4][space] Mean=1200] ] |
| ] |
| |
| Note that for large /x/ these functions are largely dependent on |
| the accuracy of the `std::sin` and `std::cos` functions. |
| |
| Comparison to GSL and __cephes is interesting: both __cephes 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, __cephes 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]. Note that |
| by way of double-checking these results, the worst performing __cephes and GSL cases |
| were recomputed using [@http://functions.wolfram.com functions.wolfram.com], |
| and the result checked against our test data: no errors in the test data were found. |
| |
| [h4 Implementation] |
| |
| The implementation is mostly about filtering off various special cases: |
| |
| When /x/ is negative, then the order /v/ 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 /x > 0/. |
| |
| When the order /v/ is negative then the reflection formulae can be used to |
| move to /v > 0/: |
| |
| [equation bessel9] |
| |
| [equation bessel10] |
| |
| Note that if the order is an integer, then these formulae reduce to: |
| |
| J[sub -n] = (-1)[super n]J[sub n] |
| |
| Y[sub -n] = (-1)[super n]Y[sub n] |
| |
| However, in general, a negative order implies that we will need to compute |
| both J and Y. |
| |
| When /x/ is large compared to the order /v/ then the asymptotic expansions |
| for large /x/ in M. Abramowitz and I.A. Stegun, |
| ['Handbook of Mathematical Functions] 9.2.19 are used |
| (these were found to be more reliable |
| than those in A&S 9.2.5). |
| |
| When the order /v/ is an integer the method first relates the result |
| to J[sub 0], J[sub 1], Y[sub 0][space] and Y[sub 1][space] using either |
| forwards or backwards recurrence (Miller's algorithm) depending upon which is stable. |
| The values for J[sub 0], J[sub 1], Y[sub 0][space] and Y[sub 1][space] are |
| calculated using the rational minimax approximations on |
| root-bracketing intervals for small ['|x|] and Hankel asymptotic |
| expansion for large ['|x|]. The coefficients are from: |
| |
| W.J. Cody, ['ALGORITHM 715: SPECFUN - A Portable FORTRAN Package of |
| Special Function Routines and Test Drivers], ACM Transactions on Mathematical |
| Software, vol 19, 22 (1993). |
| |
| and |
| |
| J.F. Hart et al, ['Computer Approximations], John Wiley & Sons, New York, 1968. |
| |
| These approximations are accurate to around 19 decimal digits: therefore |
| these methods are not used when type T has more than 64 binary digits. |
| |
| When /x/ is small, J[sub x][space] is best computed directly from the series: |
| |
| [equation bessel2] |
| |
| In the general case we compute J[sub v][space] and |
| Y[sub v][space] simultaneously. |
| |
| To get the initial values, let |
| [mu][space] = [nu] - floor([nu] + 1/2), then [mu][space] is the fractional part |
| of [nu][space] such that |
| |[mu]| <= 1/2 (we need this for convergence later). The idea is to |
| calculate J[sub [mu]](x), J[sub [mu]+1](x), Y[sub [mu]](x), Y[sub [mu]+1](x) |
| and use them to obtain J[sub [nu]](x), Y[sub [nu]](x). |
| |
| The algorithm is called Steed's method, which needs two |
| continued fractions as well as the Wronskian: |
| |
| [equation bessel8] |
| |
| [equation bessel11] |
| |
| [equation bessel12] |
| |
| See: F.S. Acton, ['Numerical Methods that Work], |
| The Mathematical Association of America, Washington, 1997. |
| |
| The continued fractions are computed using the modified Lentz's method |
| (W.J. Lentz, ['Generating Bessel functions in Mie scattering calculations |
| using continued fractions], Applied Optics, vol 15, 668 (1976)). |
| Their convergence rates depend on ['x], therefore we need |
| different strategies for large ['x] and small ['x]. |
| |
| ['x > v], CF1 needs O(['x]) iterations to converge, CF2 converges rapidly |
| |
| ['x <= v], CF1 converges rapidly, CF2 fails to converge when ['x] [^->] 0 |
| |
| When ['x] is large (['x] > 2), both continued fractions converge (CF1 |
| may be slow for really large ['x]). J[sub [mu]], J[sub [mu]+1], |
| Y[sub [mu]], Y[sub [mu]+1] can be calculated by |
| |
| [equation bessel13] |
| |
| where |
| |
| [equation bessel14] |
| |
| J[sub [nu]] and Y[sub [mu]] are then calculated using backward |
| (Miller's algorithm) and forward recurrence respectively. |
| |
| When ['x] is small (['x] <= 2), CF2 convergence may fail (but CF1 |
| works very well). The solution here is Temme's series: |
| |
| [equation bessel15] |
| |
| where |
| |
| [equation bessel16] |
| |
| g[sub k][space] and h[sub k][space] |
| are also computed by recursions (involving gamma functions), but the |
| formulas are a little complicated, readers are refered to |
| N.M. Temme, ['On the numerical evaluation of the ordinary Bessel function |
| of the second kind], Journal of Computational Physics, vol 21, 343 (1976). |
| Note Temme's series converge only for |[mu]| <= 1/2. |
| |
| As the previous case, Y[sub [nu]][space] is calculated from the forward |
| recurrence, so is Y[sub [nu]+1]. With these two |
| values and f[sub [nu]], the Wronskian yields J[sub [nu]](x) directly |
| without backward recurrence. |
| |
| [endsect] |
| |
| [/ |
| Copyright 2006 John Maddock, Paul A. Bristow and Xiaogang Zhang. |
| Distributed under the Boost Software License, Version 1.0. |
| (See accompanying file LICENSE_1_0.txt or copy at |
| http://www.boost.org/LICENSE_1_0.txt). |
| ] |