]>
Commit | Line | Data |
---|---|---|
7c673cae FG |
1 | |
2 | [section:bessel_first Bessel Functions of the First and Second Kinds] | |
3 | ||
4 | [h4 Synopsis] | |
5 | ||
6 | `#include <boost/math/special_functions/bessel.hpp>` | |
7 | ||
8 | template <class T1, class T2> | |
9 | ``__sf_result`` cyl_bessel_j(T1 v, T2 x); | |
10 | ||
11 | template <class T1, class T2, class ``__Policy``> | |
12 | ``__sf_result`` cyl_bessel_j(T1 v, T2 x, const ``__Policy``&); | |
13 | ||
14 | template <class T1, class T2> | |
15 | ``__sf_result`` cyl_neumann(T1 v, T2 x); | |
16 | ||
17 | template <class T1, class T2, class ``__Policy``> | |
18 | ``__sf_result`` cyl_neumann(T1 v, T2 x, const ``__Policy``&); | |
19 | ||
20 | ||
21 | [h4 Description] | |
22 | ||
23 | The functions __cyl_bessel_j and __cyl_neumann return the result of the | |
24 | Bessel functions of the first and second kinds respectively: | |
25 | ||
26 | cyl_bessel_j(v, x) = J[sub v](x) | |
27 | ||
28 | cyl_neumann(v, x) = Y[sub v](x) = N[sub v](x) | |
29 | ||
30 | where: | |
31 | ||
32 | [equation bessel2] | |
33 | ||
34 | [equation bessel3] | |
35 | ||
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. | |
39 | ||
40 | [optional_policy] | |
41 | ||
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 | |
45 | when `x <= 0`. | |
46 | ||
47 | The following graph illustrates the cyclic nature of J[sub v]: | |
48 | ||
49 | [graph cyl_bessel_j] | |
50 | ||
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/: | |
53 | ||
54 | [graph cyl_neumann] | |
55 | ||
56 | [h4 Testing] | |
57 | ||
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). | |
63 | ||
64 | [h4 Accuracy] | |
65 | ||
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 | |
78 | information on this. | |
79 | ||
80 | [table_cyl_bessel_j_integer_orders_] | |
81 | ||
82 | [table_cyl_bessel_j] | |
83 | ||
84 | [table_cyl_neumann_integer_orders_] | |
85 | ||
86 | [table_cyl_neumann] | |
87 | ||
88 | Note that for large /x/ these functions are largely dependent on | |
89 | the accuracy of the `std::sin` and `std::cos` functions. | |
90 | ||
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. | |
101 | ||
102 | [h4 Implementation] | |
103 | ||
104 | The implementation is mostly about filtering off various special cases: | |
105 | ||
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/. | |
109 | ||
110 | When the order /v/ is negative then the reflection formulae can be used to | |
111 | move to /v > 0/: | |
112 | ||
113 | [equation bessel9] | |
114 | ||
115 | [equation bessel10] | |
116 | ||
117 | Note that if the order is an integer, then these formulae reduce to: | |
118 | ||
119 | J[sub -n] = (-1)[super n]J[sub n] | |
120 | ||
121 | Y[sub -n] = (-1)[super n]Y[sub n] | |
122 | ||
123 | However, in general, a negative order implies that we will need to compute | |
124 | both J and Y. | |
125 | ||
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). | |
131 | ||
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: | |
139 | ||
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). | |
143 | ||
144 | and | |
145 | ||
146 | J.F. Hart et al, ['Computer Approximations], John Wiley & Sons, New York, 1968. | |
147 | ||
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. | |
150 | ||
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]): | |
157 | ||
158 | [equation bessel_y0_small_z] | |
159 | ||
160 | [equation bessel_y1_small_z] | |
161 | ||
162 | [equation bessel_y2_small_z] | |
163 | ||
164 | [equation bessel_yn_small_z] | |
165 | ||
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]): | |
170 | ||
171 | [equation bessel_yv_small_z] | |
172 | ||
173 | When /x/ is small compared to /v/, J[sub v]x[space] is best computed directly from the series: | |
174 | ||
175 | [equation bessel2] | |
176 | ||
177 | In the general case we compute J[sub v][space] and | |
178 | Y[sub v][space] simultaneously. | |
179 | ||
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). | |
186 | ||
187 | The algorithm is called Steed's method, which needs two | |
188 | continued fractions as well as the Wronskian: | |
189 | ||
190 | [equation bessel8] | |
191 | ||
192 | [equation bessel11] | |
193 | ||
194 | [equation bessel12] | |
195 | ||
196 | See: F.S. Acton, ['Numerical Methods that Work], | |
197 | The Mathematical Association of America, Washington, 1997. | |
198 | ||
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]. | |
204 | ||
205 | ['x > v], CF1 needs O(['x]) iterations to converge, CF2 converges rapidly | |
206 | ||
207 | ['x <= v], CF1 converges rapidly, CF2 fails to converge when ['x] [^->] 0 | |
208 | ||
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 | |
212 | ||
213 | [equation bessel13] | |
214 | ||
215 | where | |
216 | ||
217 | [equation bessel14] | |
218 | ||
219 | J[sub [nu]] and Y[sub [mu]] are then calculated using backward | |
220 | (Miller's algorithm) and forward recurrence respectively. | |
221 | ||
222 | When ['x] is small (['x] <= 2), CF2 convergence may fail (but CF1 | |
223 | works very well). The solution here is Temme's series: | |
224 | ||
225 | [equation bessel15] | |
226 | ||
227 | where | |
228 | ||
229 | [equation bessel16] | |
230 | ||
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. | |
237 | ||
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. | |
242 | ||
243 | [endsect] | |
244 | ||
245 | [section:bessel_root Finding Zeros of Bessel Functions of the First and Second Kinds] | |
246 | ||
247 | [h4 Synopsis] | |
248 | ||
249 | `#include <boost/math/special_functions/bessel.hpp>` | |
250 | ||
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. | |
254 | ||
255 | The signature of the single value functions are: | |
256 | ||
257 | template <class T> | |
258 | T cyl_bessel_j_zero( | |
259 | T v, // Floating-point value for Jv. | |
260 | int m); // 1-based index of zero. | |
261 | ||
262 | template <class T> | |
263 | T cyl_neumann_zero( | |
264 | T v, // Floating-point value for Jv. | |
265 | int m); // 1-based index of zero. | |
266 | ||
267 | and for multiple zeros: | |
268 | ||
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. | |
275 | ||
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. | |
282 | ||
283 | There are also versions which allow control of the __policy_section for error handling and precision. | |
284 | ||
285 | template <class T> | |
286 | T cyl_bessel_j_zero( | |
287 | T v, // Floating-point value for Jv. | |
288 | int m, // 1-based index of zero. | |
289 | const Policy&); // Policy to use. | |
290 | ||
291 | template <class T> | |
292 | T cyl_neumann_zero( | |
293 | T v, // Floating-point value for Jv. | |
294 | int m, // 1-based index of zero. | |
295 | const Policy&); // Policy to use. | |
296 | ||
297 | ||
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. | |
305 | ||
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. | |
313 | ||
314 | [h4 Description] | |
315 | ||
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 | |
319 | ||
320 | [emquad] ['J[sub [nu]](j[sub [nu], m]) = 0] | |
321 | ||
322 | [emquad] ['Y[sub [nu]](y[sub [nu], m]) = 0] | |
323 | ||
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]]. | |
328 | ||
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`. | |
332 | ||
333 | In each case the index or rank of the zero | |
334 | returned is 1-based, which is to say: | |
335 | ||
336 | cyl_bessel_j_zero(v, 1); | |
337 | ||
338 | returns the first zero of Bessel J. | |
339 | ||
340 | Passing an `start_index <= 0` results in a `std::domain_error` being raised. | |
341 | ||
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. | |
347 | ||
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. | |
350 | ||
351 | [graph bessel_j_zeros] | |
352 | ||
353 | [graph neumann_y_zeros] | |
354 | ||
355 | [h4 Examples of finding Bessel and Neumann zeros] | |
356 | ||
357 | [import ../../example/bessel_zeros_example_1.cpp] | |
358 | ||
359 | [bessel_zeros_example_1] | |
360 | [bessel_zeros_example_2] | |
361 | ||
362 | [import ../../example/bessel_zeros_interator_example.cpp] | |
363 | ||
364 | [bessel_zeros_iterator_example_1] | |
365 | [bessel_zeros_iterator_example_2] | |
366 | ||
367 | [import ../../example/neumann_zeros_example_1.cpp] | |
368 | ||
369 | [neumann_zeros_example_1] | |
370 | [neumann_zeros_example_2] | |
371 | ||
372 | [import ../../example/bessel_errors_example.cpp] | |
373 | ||
374 | [bessel_errors_example_1] | |
375 | [bessel_errors_example_2] | |
376 | ||
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]. | |
382 | ||
383 | [h3 Implementation] | |
384 | ||
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. | |
387 | ||
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. | |
392 | ||
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. | |
396 | In particular, | |
397 | ||
398 | [emquad] ['d/[sub dx] ['J[sub [nu]](x)] = ['J[sub [nu]-1](x)] - [nu] J[sub [nu]](x)] / x | |
399 | ||
400 | [emquad] ['d/[sub dx] ['Y[sub [nu]](x)] = ['Y[sub [nu]-1](x)] - [nu] Y[sub [nu]](x)] / x | |
401 | ||
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. | |
405 | ||
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]. | |
411 | ||
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 | |
414 | and ['n [ne] 0]. | |
415 | ||
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()`. | |
420 | ||
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. | |
430 | ||
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. | |
435 | ||
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. | |
438 | ||
439 | In particular, | |
440 | ||
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]] | |
442 | ||
443 | and | |
444 | ||
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]] | |
446 | ||
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. | |
450 | In particular, | |
451 | ||
452 | [emquad] ['j[sub [nu],m], y[sub [nu],m] [cong] [beta] - ([mu]-1) / 8[beta]] | |
453 | ||
454 | [emquad] [emquad] [emquad] ['- 4([mu]-1)(7[mu] - 31) / 3(8[beta])[super 3]] | |
455 | ||
456 | [emquad] [emquad] [emquad] ['-32([mu]-1)(83[mu][sup2] - 982[mu] + 3779) / 15(8[beta])[super 5]] | |
457 | ||
458 | [emquad] [emquad] [emquad] ['-64([mu]-1)(6949[mu][super 3] - 153855[mu][sup2] + 1585743[mu]- 6277237) / 105(8a)[super 7]] | |
459 | ||
460 | [emquad] [emquad] [emquad] ['- [ellipsis]] [emquad] [emquad] (5) | |
461 | ||
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]]]. | |
465 | ||
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). | |
475 | ||
476 | In summary, | |
477 | ||
478 | [emquad] ['j[sub [nu], m] [sim] [nu]x(-[zeta]) + f[sub 1](-[zeta]/[nu])] | |
479 | ||
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. | |
482 | ||
483 | Here ['x = x(-[zeta])] is the inverse of the function | |
484 | ||
485 | [emquad] ['[frac23](-[zeta])[super 3/2] = [radic](x[sup2] - 1) - cos[supminus][sup1](1/x)] [emquad] [emquad] (7) | |
486 | ||
487 | Furthermore, | |
488 | ||
489 | [emquad] ['f[sub 1](-[zeta]) = [frac12]x(-[zeta]) {h(-[zeta])}[sup2] [sdot] b[sub 0](-[zeta])] | |
490 | ||
491 | where | |
492 | ||
493 | [emquad] ['h(-[zeta]) = {4(-[zeta]) / (x[sup2] - 1)}[super 4]] | |
494 | ||
495 | and | |
496 | ||
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])]}] | |
498 | ||
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 | |
502 | ||
503 | [emquad] ['[frac23](-[zeta])[super 3/2] [approx] x + 1/2x - [pi]/2] | |
504 | ||
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 | |
509 | in Eq. 7. | |
510 | ||
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, | |
518 | ||
519 | [emquad] ['j[sub n[nu],m]] < ['j[sub -[nu],m]] < ['j[sub n[nu],m+1]] | |
520 | ||
521 | where ['m > 1] and ['n[sub [nu]]] represents the integral | |
522 | floor of the absolute value of ['|-[nu]|]. | |
523 | ||
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. | |
528 | ||
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]]. | |
534 | ||
535 | [h3 Testing] | |
536 | ||
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. | |
539 | ||
540 | [endsect] [/section:bessel Finding Zeros of Bessel Functions of the First and Second Kinds] | |
541 | ||
542 | [/ | |
543 | Copyright 2006, 2013 John Maddock, Paul A. Bristow, Xiaogang Zhang and Christopher Kormanyos. | |
544 | ||
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). | |
548 | ] |