libs/math/doc/sf_and_dist/html/math_toolkit/special/bessel/bessel.html

   1 <html>
   2 <head>
   3 <meta http-equiv="Content-Type" content="text/html; charset=US-ASCII">
   4 <title>Bessel Functions of the First and Second Kinds</title>
   5 <link rel="stylesheet" href="../../../../../../../../doc/src/boostbook.css" type="text/css">
   6 <meta name="generator" content="DocBook XSL Stylesheets V1.76.1">
   7 <link rel="home" href="../../../index.html" title="Math Toolkit">
   8 <link rel="up" href="../bessel.html" title="Bessel Functions">
   9 <link rel="prev" href="bessel_over.html" title="Bessel Function Overview">
  10 <link rel="next" href="mbessel.html" title="Modified Bessel Functions of the First and Second Kinds">
  11 </head>
  12 <body bgcolor="white" text="black" link="#0000FF" vlink="#840084" alink="#0000FF">
  13 <table cellpadding="2" width="100%"><tr>
  14 <td valign="top"><img alt="Boost C++ Libraries" width="277" height="86" src="../../../../../../../../boost.png"></td>
  15 <td align="center"><a href="../../../../../../../../index.html">Home</a></td>
  16 <td align="center"><a href="../../../../../../../../libs/libraries.htm">Libraries</a></td>
  17 <td align="center"><a href="http://www.boost.org/users/people.html">People</a></td>
  18 <td align="center"><a href="http://www.boost.org/users/faq.html">FAQ</a></td>
  19 <td align="center"><a href="../../../../../../../../more/index.htm">More</a></td>
  20 </tr></table>
  21 <hr>
  22 <div class="spirit-nav">
  23 <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="mbessel.html"><img src="../../../../../../../../doc/src/images/next.png" alt="Next"></a>
  24 </div>
  25 <div class="section math_toolkit_special_bessel_bessel">
  26 <div class="titlepage"><div><div><h4 class="title">
  27 <a name="math_toolkit.special.bessel.bessel"></a><a class="link" href="bessel.html" title="Bessel Functions of the First and Second Kinds">Bessel Functions
  28         of the First and Second Kinds</a>
  29 </h4></div></div></div>
  30 <h5>
  31 <a name="math_toolkit.special.bessel.bessel.h0"></a>
  32           <span><a name="math_toolkit.special.bessel.bessel.synopsis"></a></span><a class="link" href="bessel.html#math_toolkit.special.bessel.bessel.synopsis">Synopsis</a>
  33         </h5>
  34 <pre class="programlisting"><span class="keyword">template</span> <span class="special">&lt;</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">&gt;</span>
  35 <a class="link" href="../../main_overview/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>
  36
  37 <span class="keyword">template</span> <span class="special">&lt;</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="Policies">Policy</a><span class="special">&gt;</span>
  38 <a class="link" href="../../main_overview/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="Policies">Policy</a><span class="special">&amp;);</span>
  39
  40 <span class="keyword">template</span> <span class="special">&lt;</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">&gt;</span>
  41 <a class="link" href="../../main_overview/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>
  42
  43 <span class="keyword">template</span> <span class="special">&lt;</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="Policies">Policy</a><span class="special">&gt;</span>
  44 <a class="link" href="../../main_overview/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="Policies">Policy</a><span class="special">&amp;);</span>
  45 </pre>
  46 <h5>
  47 <a name="math_toolkit.special.bessel.bessel.h1"></a>
  48           <span><a name="math_toolkit.special.bessel.bessel.description"></a></span><a class="link" href="bessel.html#math_toolkit.special.bessel.bessel.description">Description</a>
  49         </h5>
  50 <p>
  51           The functions <a class="link" href="bessel.html" title="Bessel Functions of the First and Second Kinds">cyl_bessel_j</a>
  52           and <a class="link" href="bessel.html" title="Bessel Functions of the First and Second Kinds">cyl_neumann</a>
  53           return the result of the Bessel functions of the first and second kinds
  54           respectively:
  55         </p>
  56 <p>
  57           cyl_bessel_j(v, x) = J<sub>v</sub>(x)
  58         </p>
  59 <p>
  60           cyl_neumann(v, x) = Y<sub>v</sub>(x) = N<sub>v</sub>(x)
  61         </p>
  62 <p>
  63           where:
  64         </p>
  65 <p>
  66           <span class="inlinemediaobject"><img src="../../../../equations/bessel2.png"></span>
  67         </p>
  68 <p>
  69           <span class="inlinemediaobject"><img src="../../../../equations/bessel3.png"></span>
  70         </p>
  71 <p>
  72           The return type of these functions is computed using the <a class="link" href="../../main_overview/result_type.html" title="Calculation of the Type of the Result"><span class="emphasis"><em>result
  73           type calculation rules</em></span></a> when T1 and T2 are different types.
  74           The functions are also optimised for the relatively common case that T1
  75           is an integer.
  76         </p>
  77 <p>
  78           The final <a class="link" href="../../policy.html" title="Policies">Policy</a> argument is
  79           optional and can be used to control the behaviour of the function: how
  80           it handles errors, what level of precision to use etc. Refer to the <a class="link" href="../../policy.html" title="Policies">policy documentation for more details</a>.
  81         </p>
  82 <p>
  83           The functions return the result of <a class="link" href="../../main_overview/error_handling.html#domain_error">domain_error</a>
  84           whenever the result is undefined or complex. For <a class="link" href="bessel.html" title="Bessel Functions of the First and Second Kinds">cyl_bessel_j</a>
  85           this occurs when <code class="computeroutput"><span class="identifier">x</span> <span class="special">&lt;</span>
  86           <span class="number">0</span></code> and v is not an integer, or when
  87           <code class="computeroutput"><span class="identifier">x</span> <span class="special">==</span>
  88           <span class="number">0</span></code> and <code class="computeroutput"><span class="identifier">v</span>
  89           <span class="special">!=</span> <span class="number">0</span></code>.
  90           For <a class="link" href="bessel.html" title="Bessel Functions of the First and Second Kinds">cyl_neumann</a>
  91           this occurs when <code class="computeroutput"><span class="identifier">x</span> <span class="special">&lt;=</span>
  92           <span class="number">0</span></code>.
  93         </p>
  94 <p>
  95           The following graph illustrates the cyclic nature of J<sub>v</sub>:
  96         </p>
  97 <p>
  98           <span class="inlinemediaobject"><img src="../../../../graphs/cyl_bessel_j.png" align="middle"></span>
  99         </p>
 100 <p>
 101           The following graph shows the behaviour of Y<sub>v</sub>: this is also cyclic for
 102           large <span class="emphasis"><em>x</em></span>, but tends to -&#8734; &#160; for small <span class="emphasis"><em>x</em></span>:
 103         </p>
 104 <p>
 105           <span class="inlinemediaobject"><img src="../../../../graphs/cyl_neumann.png" align="middle"></span>
 106         </p>
 107 <h5>
 108 <a name="math_toolkit.special.bessel.bessel.h2"></a>
 109           <span><a name="math_toolkit.special.bessel.bessel.testing"></a></span><a class="link" href="bessel.html#math_toolkit.special.bessel.bessel.testing">Testing</a>
 110         </h5>
 111 <p>
 112           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
 113           much larger set of tests computed using a simplified version of this implementation
 114           (with all the special case handling removed).
 115         </p>
 116 <h5>
 117 <a name="math_toolkit.special.bessel.bessel.h3"></a>
 118           <span><a name="math_toolkit.special.bessel.bessel.accuracy"></a></span><a class="link" href="bessel.html#math_toolkit.special.bessel.bessel.accuracy">Accuracy</a>
 119         </h5>
 120 <p>
 121           The following tables show how the accuracy of these functions varies on
 122           various platforms, along with comparisons to the <a href="http://www.gnu.org/software/gsl/" target="_top">GSL-1.9</a>
 123           and <a href="http://www.netlib.org/cephes/" target="_top">Cephes</a> libraries.
 124           Note that the cyclic nature of these functions means that they have an
 125           infinite number of irrational roots: in general these functions have arbitrarily
 126           large <span class="emphasis"><em>relative</em></span> errors when the arguments are sufficiently
 127           close to a root. Of course the absolute error in such cases is always small.
 128           Note that only results for the widest floating-point type on the system
 129           are given as narrower types have <a class="link" href="../../backgrounders/relative_error.html#zero_error">effectively
 130           zero error</a>. All values are relative errors in units of epsilon.
 131         </p>
 132 <div class="table">
 133 <a name="math_toolkit.special.bessel.bessel.errors_rates_in_cyl_bessel_j"></a><p class="title"><b>Table&#160;39.&#160;Errors Rates in cyl_bessel_j</b></p>
 134 <div class="table-contents"><table class="table" summary="Errors Rates in cyl_bessel_j">
 135 <colgroup>
 136 <col>
 137 <col>
 138 <col>
 139 <col>
 140 <col>
 141 </colgroup>
 142 <thead><tr>
 143 <th>
 144                   <p>
 145                     Significand Size
 146                   </p>
 147                 </th>
 148 <th>
 149                   <p>
 150                     Platform and Compiler
 151                   </p>
 152                 </th>
 153 <th>
 154                   <p>
 155                     J<sub>0</sub> &#160; and J<sub>1</sub>
 156                   </p>
 157                 </th>
 158 <th>
 159                   <p>
 160                     J<sub>v</sub>
 161                   </p>
 162                 </th>
 163 <th>
 164                   <p>
 165                     J<sub>v</sub> &#160; (large values of x &gt; 1000)
 166                   </p>
 167                 </th>
 168 </tr></thead>
 169 <tbody>
 170 <tr>
 171 <td>
 172                   <p>
 173                     53
 174                   </p>
 175                 </td>
 176 <td>
 177                   <p>
 178                     Win32 / Visual C++ 8.0
 179                   </p>
 180                 </td>
 181 <td>
 182                   <p>
 183                     Peak=2.5 Mean=1.1
 184                   </p>
 185                   <p>
 186                     GSL Peak=6.6
 187                   </p>
 188                   <p>
 189                     <a href="http://www.netlib.org/cephes/" target="_top">Cephes</a> Peak=2.5
 190                     Mean=1.1
 191                   </p>
 192                 </td>
 193 <td>
 194                   <p>
 195                     Peak=11 Mean=2.2
 196                   </p>
 197                   <p>
 198                     GSL Peak=11
 199                   </p>
 200                   <p>
 201                     <a href="http://www.netlib.org/cephes/" target="_top">Cephes</a> Peak=17
 202                     Mean=2.5
 203                   </p>
 204                 </td>
 205 <td>
 206                   <p>
 207                     Peak=413 Mean=110
 208                   </p>
 209                   <p>
 210                     GSL Peak=6x10<sup>11</sup>
 211                   </p>
 212                   <p>
 213                     <a href="http://www.netlib.org/cephes/" target="_top">Cephes</a> Peak=2x10<sup>5</sup>
 214                   </p>
 215                 </td>
 216 </tr>
 217 <tr>
 218 <td>
 219                   <p>
 220                     64
 221                   </p>
 222                 </td>
 223 <td>
 224                   <p>
 225                     Red Hat Linux IA64 / G++ 3.4
 226                   </p>
 227                 </td>
 228 <td>
 229                   <p>
 230                     Peak=7 Mean=3
 231                   </p>
 232                 </td>
 233 <td>
 234                   <p>
 235                     Peak=117 Mean=10
 236                   </p>
 237                 </td>
 238 <td>
 239                   <p>
 240                     Peak=2x10<sup>4</sup> &#160; Mean=6x10<sup>3</sup>
 241                   </p>
 242                 </td>
 243 </tr>
 244 <tr>
 245 <td>
 246                   <p>
 247                     64
 248                   </p>
 249                 </td>
 250 <td>
 251                   <p>
 252                     SUSE Linux AMD64 / G++ 4.1
 253                   </p>
 254                 </td>
 255 <td>
 256                   <p>
 257                     Peak=7 Mean=3
 258                   </p>
 259                 </td>
 260 <td>
 261                   <p>
 262                     Peak=400 Mean=40
 263                   </p>
 264                 </td>
 265 <td>
 266                   <p>
 267                     Peak=2x10<sup>4</sup> &#160; Mean=1x10<sup>4</sup>
 268                   </p>
 269                 </td>
 270 </tr>
 271 <tr>
 272 <td>
 273                   <p>
 274                     113
 275                   </p>
 276                 </td>
 277 <td>
 278                   <p>
 279                     HP-UX / HP aCC 6
 280                   </p>
 281                 </td>
 282 <td>
 283                   <p>
 284                     Peak=14 Mean=6
 285                   </p>
 286                 </td>
 287 <td>
 288                   <p>
 289                     Peak=29 Mean=3
 290                   </p>
 291                 </td>
 292 <td>
 293                   <p>
 294                     Peak=2700 Mean=450
 295                   </p>
 296                 </td>
 297 </tr>
 298 </tbody>
 299 </table></div>
 300 </div>
 301 <br class="table-break"><div class="table">
 302 <a name="math_toolkit.special.bessel.bessel.errors_rates_in_cyl_neumann"></a><p class="title"><b>Table&#160;40.&#160;Errors Rates in cyl_neumann</b></p>
 303 <div class="table-contents"><table class="table" summary="Errors Rates in cyl_neumann">
 304 <colgroup>
 305 <col>
 306 <col>
 307 <col>
 308 <col>
 309 <col>
 310 </colgroup>
 311 <thead><tr>
 312 <th>
 313                   <p>
 314                     Significand Size
 315                   </p>
 316                 </th>
 317 <th>
 318                   <p>
 319                     Platform and Compiler
 320                   </p>
 321                 </th>
 322 <th>
 323                   <p>
 324                     J<sub>0</sub> &#160; and J<sub>1</sub>
 325                   </p>
 326                 </th>
 327 <th>
 328                   <p>
 329                     J<sub>n</sub> (integer orders)
 330                   </p>
 331                 </th>
 332 <th>
 333                   <p>
 334                     J<sub>v</sub> (fractional orders)
 335                   </p>
 336                 </th>
 337 </tr></thead>
 338 <tbody>
 339 <tr>
 340 <td>
 341                   <p>
 342                     53
 343                   </p>
 344                 </td>
 345 <td>
 346                   <p>
 347                     Win32 / Visual C++ 8.0
 348                   </p>
 349                 </td>
 350 <td>
 351                   <p>
 352                     Peak=330 Mean=54
 353                   </p>
 354                   <p>
 355                     GSL Peak=34 Mean=9
 356                   </p>
 357                   <p>
 358                     <a href="http://www.netlib.org/cephes/" target="_top">Cephes</a> Peak=330
 359                     Mean=54
 360                   </p>
 361                 </td>
 362 <td>
 363                   <p>
 364                     Peak=923 Mean=83
 365                   </p>
 366                   <p>
 367                     GSL Peak=500 Mean=54
 368                   </p>
 369                   <p>
 370                     <a href="http://www.netlib.org/cephes/" target="_top">Cephes</a> Peak=923
 371                     Mean=83
 372                   </p>
 373                 </td>
 374 <td>
 375                   <p>
 376                     Peak=561 Mean=36
 377                   </p>
 378                   <p>
 379                     GSL Peak=1.4x10<sup>6</sup> &#160; Mean=7x10<sup>4</sup> &#160;
 380                   </p>
 381                   <p>
 382                     <a href="http://www.netlib.org/cephes/" target="_top">Cephes</a> Peak=+INF
 383                   </p>
 384                 </td>
 385 </tr>
 386 <tr>
 387 <td>
 388                   <p>
 389                     64
 390                   </p>
 391                 </td>
 392 <td>
 393                   <p>
 394                     Red Hat Linux IA64 / G++ 3.4
 395                   </p>
 396                 </td>
 397 <td>
 398                   <p>
 399                     Peak=470 Mean=56
 400                   </p>
 401                 </td>
 402 <td>
 403                   <p>
 404                     Peak=843 Mean=51
 405                   </p>
 406                 </td>
 407 <td>
 408                   <p>
 409                     Peak=741 Mean=51
 410                   </p>
 411                 </td>
 412 </tr>
 413 <tr>
 414 <td>
 415                   <p>
 416                     64
 417                   </p>
 418                 </td>
 419 <td>
 420                   <p>
 421                     SUSE Linux AMD64 / G++ 4.1
 422                   </p>
 423                 </td>
 424 <td>
 425                   <p>
 426                     Peak=1300 Mean=424
 427                   </p>
 428                 </td>
 429 <td>
 430                   <p>
 431                     Peak=2x10<sup>4</sup> &#160; Mean=8x10<sup>3</sup>
 432                   </p>
 433                 </td>
 434 <td>
 435                   <p>
 436                     Peak=1x10<sup>5</sup> &#160; Mean=6x10<sup>3</sup>
 437                   </p>
 438                 </td>
 439 </tr>
 440 <tr>
 441 <td>
 442                   <p>
 443                     113
 444                   </p>
 445                 </td>
 446 <td>
 447                   <p>
 448                     HP-UX / HP aCC 6
 449                   </p>
 450                 </td>
 451 <td>
 452                   <p>
 453                     Peak=180 Mean=63
 454                   </p>
 455                 </td>
 456 <td>
 457                   <p>
 458                     Peak=340 Mean=150
 459                   </p>
 460                 </td>
 461 <td>
 462                   <p>
 463                     Peak=2x10<sup>4</sup> &#160; Mean=1200
 464                   </p>
 465                 </td>
 466 </tr>
 467 </tbody>
 468 </table></div>
 469 </div>
 470 <br class="table-break"><p>
 471           Note that for large <span class="emphasis"><em>x</em></span> these functions are largely
 472           dependent on the accuracy of the <code class="computeroutput"><span class="identifier">std</span><span class="special">::</span><span class="identifier">sin</span></code>
 473           and <code class="computeroutput"><span class="identifier">std</span><span class="special">::</span><span class="identifier">cos</span></code> functions.
 474         </p>
 475 <p>
 476           Comparison to GSL and <a href="http://www.netlib.org/cephes/" target="_top">Cephes</a>
 477           is interesting: both <a href="http://www.netlib.org/cephes/" target="_top">Cephes</a>
 478           and this library optimise the integer order case - leading to identical
 479           results - simply using the general case is for the most part slightly more
 480           accurate though, as noted by the better accuracy of GSL in the integer
 481           argument cases. This implementation tends to perform much better when the
 482           arguments become large, <a href="http://www.netlib.org/cephes/" target="_top">Cephes</a>
 483           in particular produces some remarkably inaccurate results with some of
 484           the test data (no significant figures correct), and even GSL performs badly
 485           with some inputs to J<sub>v</sub>. Note that by way of double-checking these results,
 486           the worst performing <a href="http://www.netlib.org/cephes/" target="_top">Cephes</a>
 487           and GSL cases were recomputed using <a href="http://functions.wolfram.com" target="_top">functions.wolfram.com</a>,
 488           and the result checked against our test data: no errors in the test data
 489           were found.
 490         </p>
 491 <h5>
 492 <a name="math_toolkit.special.bessel.bessel.h4"></a>
 493           <span><a name="math_toolkit.special.bessel.bessel.implementation"></a></span><a class="link" href="bessel.html#math_toolkit.special.bessel.bessel.implementation">Implementation</a>
 494         </h5>
 495 <p>
 496           The implementation is mostly about filtering off various special cases:
 497         </p>
 498 <p>
 499           When <span class="emphasis"><em>x</em></span> is negative, then the order <span class="emphasis"><em>v</em></span>
 500           must be an integer or the result is a domain error. If the order is an
 501           integer then the function is odd for odd orders and even for even orders,
 502           so we reflect to <span class="emphasis"><em>x &gt; 0</em></span>.
 503         </p>
 504 <p>
 505           When the order <span class="emphasis"><em>v</em></span> is negative then the reflection formulae
 506           can be used to move to <span class="emphasis"><em>v &gt; 0</em></span>:
 507         </p>
 508 <p>
 509           <span class="inlinemediaobject"><img src="../../../../equations/bessel9.png"></span>
 510         </p>
 511 <p>
 512           <span class="inlinemediaobject"><img src="../../../../equations/bessel10.png"></span>
 513         </p>
 514 <p>
 515           Note that if the order is an integer, then these formulae reduce to:
 516         </p>
 517 <p>
 518           J<sub>-n</sub> = (-1)<sup>n</sup>J<sub>n</sub>
 519         </p>
 520 <p>
 521           Y<sub>-n</sub> = (-1)<sup>n</sup>Y<sub>n</sub>
 522         </p>
 523 <p>
 524           However, in general, a negative order implies that we will need to compute
 525           both J and Y.
 526         </p>
 527 <p>
 528           When <span class="emphasis"><em>x</em></span> is large compared to the order <span class="emphasis"><em>v</em></span>
 529           then the asymptotic expansions for large <span class="emphasis"><em>x</em></span> in M. Abramowitz
 530           and I.A. Stegun, <span class="emphasis"><em>Handbook of Mathematical Functions</em></span>
 531           9.2.19 are used (these were found to be more reliable than those in A&amp;S
 532           9.2.5).
 533         </p>
 534 <p>
 535           When the order <span class="emphasis"><em>v</em></span> is an integer the method first relates
 536           the result to J<sub>0</sub>, J<sub>1</sub>, Y<sub>0</sub> &#160; and Y<sub>1</sub> &#160; using either forwards or backwards recurrence
 537           (Miller's algorithm) depending upon which is stable. The values for J<sub>0</sub>,
 538           J<sub>1</sub>, Y<sub>0</sub> &#160; and Y<sub>1</sub> &#160; are calculated using the rational minimax approximations on
 539           root-bracketing intervals for small <span class="emphasis"><em>|x|</em></span> and Hankel
 540           asymptotic expansion for large <span class="emphasis"><em>|x|</em></span>. The coefficients
 541           are from:
 542         </p>
 543 <p>
 544           W.J. Cody, <span class="emphasis"><em>ALGORITHM 715: SPECFUN - A Portable FORTRAN Package
 545           of Special Function Routines and Test Drivers</em></span>, ACM Transactions
 546           on Mathematical Software, vol 19, 22 (1993).
 547         </p>
 548 <p>
 549           and
 550         </p>
 551 <p>
 552           J.F. Hart et al, <span class="emphasis"><em>Computer Approximations</em></span>, John Wiley
 553           &amp; Sons, New York, 1968.
 554         </p>
 555 <p>
 556           These approximations are accurate to around 19 decimal digits: therefore
 557           these methods are not used when type T has more than 64 binary digits.
 558         </p>
 559 <p>
 560           When <span class="emphasis"><em>x</em></span> is smaller than machine epsilon then the following
 561           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>,
 562           <a href="http://functions.wolfram.com/03.03.06.0038.01" target="_top">http://functions.wolfram.com/03.03.06.0038.01</a>,
 563           <a href="http://functions.wolfram.com/03.03.06.0039.01" target="_top">http://functions.wolfram.com/03.03.06.0039.01</a>
 564           and <a href="http://functions.wolfram.com/03.03.06.0040.01" target="_top">http://functions.wolfram.com/03.03.06.0040.01</a>):
 565         </p>
 566 <p>
 567           <span class="inlinemediaobject"><img src="../../../../equations/bessel_y0_small_z.png"></span>
 568         </p>
 569 <p>
 570           <span class="inlinemediaobject"><img src="../../../../equations/bessel_y1_small_z.png"></span>
 571         </p>
 572 <p>
 573           <span class="inlinemediaobject"><img src="../../../../equations/bessel_y2_small_z.png"></span>
 574         </p>
 575 <p>
 576           <span class="inlinemediaobject"><img src="../../../../equations/bessel_yn_small_z.png"></span>
 577         </p>
 578 <p>
 579           When <span class="emphasis"><em>x</em></span> is small compared to <span class="emphasis"><em>v</em></span>
 580           and <span class="emphasis"><em>v</em></span> is not an integer, then the following series
 581           approximation can be used for Y<sub>v</sub>(x), this is also an area where other approximations
 582           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>):
 583         </p>
 584 <p>
 585           <span class="inlinemediaobject"><img src="../../../../equations/bessel_yv_small_z.png"></span>
 586         </p>
 587 <p>
 588           When <span class="emphasis"><em>x</em></span> is small compared to <span class="emphasis"><em>v</em></span>,
 589           J<sub>v</sub>x &#160; is best computed directly from the series:
 590         </p>
 591 <p>
 592           <span class="inlinemediaobject"><img src="../../../../equations/bessel2.png"></span>
 593         </p>
 594 <p>
 595           In the general case we compute J<sub>v</sub> &#160; and Y<sub>v</sub> &#160; simultaneously.
 596         </p>
 597 <p>
 598           To get the initial values, let &#956; &#160; = &#957; - floor(&#957; + 1/2), then &#956; &#160; is the fractional
 599           part of &#957; &#160; such that |&#956;| &lt;= 1/2 (we need this for convergence later). The
 600           idea is to calculate J<sub>&#956;</sub>(x), J<sub>&#956;+1</sub>(x), Y<sub>&#956;</sub>(x), Y<sub>&#956;+1</sub>(x) and use them to obtain
 601           J<sub>&#957;</sub>(x), Y<sub>&#957;</sub>(x).
 602         </p>
 603 <p>
 604           The algorithm is called Steed's method, which needs two continued fractions
 605           as well as the Wronskian:
 606         </p>
 607 <p>
 608           <span class="inlinemediaobject"><img src="../../../../equations/bessel8.png"></span>
 609         </p>
 610 <p>
 611           <span class="inlinemediaobject"><img src="../../../../equations/bessel11.png"></span>
 612         </p>
 613 <p>
 614           <span class="inlinemediaobject"><img src="../../../../equations/bessel12.png"></span>
 615         </p>
 616 <p>
 617           See: F.S. Acton, <span class="emphasis"><em>Numerical Methods that Work</em></span>, The
 618           Mathematical Association of America, Washington, 1997.
 619         </p>
 620 <p>
 621           The continued fractions are computed using the modified Lentz's method
 622           (W.J. Lentz, <span class="emphasis"><em>Generating Bessel functions in Mie scattering calculations
 623           using continued fractions</em></span>, Applied Optics, vol 15, 668 (1976)).
 624           Their convergence rates depend on <span class="emphasis"><em>x</em></span>, therefore we
 625           need different strategies for large <span class="emphasis"><em>x</em></span> and small <span class="emphasis"><em>x</em></span>.
 626         </p>
 627 <p>
 628           <span class="emphasis"><em>x &gt; v</em></span>, CF1 needs O(<span class="emphasis"><em>x</em></span>) iterations
 629           to converge, CF2 converges rapidly
 630         </p>
 631 <p>
 632           <span class="emphasis"><em>x &lt;= v</em></span>, CF1 converges rapidly, CF2 fails to converge
 633           when <span class="emphasis"><em>x</em></span> <code class="literal">-&gt;</code> 0
 634         </p>
 635 <p>
 636           When <span class="emphasis"><em>x</em></span> is large (<span class="emphasis"><em>x</em></span> &gt; 2), both
 637           continued fractions converge (CF1 may be slow for really large <span class="emphasis"><em>x</em></span>).
 638           J<sub>&#956;</sub>, J<sub>&#956;+1</sub>, Y<sub>&#956;</sub>, Y<sub>&#956;+1</sub> can be calculated by
 639         </p>
 640 <p>
 641           <span class="inlinemediaobject"><img src="../../../../equations/bessel13.png"></span>
 642         </p>
 643 <p>
 644           where
 645         </p>
 646 <p>
 647           <span class="inlinemediaobject"><img src="../../../../equations/bessel14.png"></span>
 648         </p>
 649 <p>
 650           J<sub>&#957;</sub> and Y<sub>&#956;</sub> are then calculated using backward (Miller's algorithm) and forward
 651           recurrence respectively.
 652         </p>
 653 <p>
 654           When <span class="emphasis"><em>x</em></span> is small (<span class="emphasis"><em>x</em></span> &lt;= 2),
 655           CF2 convergence may fail (but CF1 works very well). The solution here is
 656           Temme's series:
 657         </p>
 658 <p>
 659           <span class="inlinemediaobject"><img src="../../../../equations/bessel15.png"></span>
 660         </p>
 661 <p>
 662           where
 663         </p>
 664 <p>
 665           <span class="inlinemediaobject"><img src="../../../../equations/bessel16.png"></span>
 666         </p>
 667 <p>
 668           g<sub>k</sub> &#160; and h<sub>k</sub> &#160;
 669 are also computed by recursions (involving gamma functions), but
 670           the formulas are a little complicated, readers are refered to N.M. Temme,
 671           <span class="emphasis"><em>On the numerical evaluation of the ordinary Bessel function of
 672           the second kind</em></span>, Journal of Computational Physics, vol 21, 343
 673           (1976). Note Temme's series converge only for |&#956;| &lt;= 1/2.
 674         </p>
 675 <p>
 676           As the previous case, Y<sub>&#957;</sub> &#160; is calculated from the forward recurrence, so is
 677           Y<sub>&#957;+1</sub>. With these two values and f<sub>&#957;</sub>, the Wronskian yields J<sub>&#957;</sub>(x) directly without
 678           backward recurrence.
 679         </p>
 680 </div>
 681 <table xmlns:rev="http://www.cs.rpi.edu/~gregod/boost/tools/doc/revision" width="100%"><tr>
 682 <td align="left"></td>
 683 <td align="right"><div class="copyright-footer">Copyright &#169; 2006-2010 John Maddock, Paul A. Bristow, Hubert Holin, Xiaogang Zhang, Bruno
 684       Lalande, Johan R&#229;de, Gautam Sewani, Thijs van den Berg and Benjamin Sobotta<p>
 685         Distributed under the Boost Software License, Version 1.0. (See accompanying
 686         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>)
 687       </p>
 688 </div></td>
 689 </tr></table>
 690 <hr>
 691 <div class="spirit-nav">
 692 <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="mbessel.html"><img src="../../../../../../../../doc/src/images/next.png" alt="Next"></a>
 693 </div>
 694 </body>
 695 </html>