1<html> 2<head> 3<meta http-equiv="Content-Type" content="text/html; charset=UTF-8"> 4<title>Modified Bessel Functions of the First and Second Kinds</title> 5<link rel="stylesheet" href="../../math.css" type="text/css"> 6<meta name="generator" content="DocBook XSL Stylesheets V1.79.1"> 7<link rel="home" href="../../index.html" title="Math Toolkit 2.12.0"> 8<link rel="up" href="../bessel.html" title="Bessel Functions"> 9<link rel="prev" href="bessel_root.html" title="Finding Zeros of Bessel Functions of the First and Second Kinds"> 10<link rel="next" href="sph_bessel.html" title="Spherical 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_root.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="sph_bessel.html"><img src="../../../../../../doc/src/images/next.png" alt="Next"></a> 24</div> 25<div class="section"> 26<div class="titlepage"><div><div><h3 class="title"> 27<a name="math_toolkit.bessel.mbessel"></a><a class="link" href="mbessel.html" title="Modified Bessel Functions of the First and Second Kinds">Modified Bessel Functions 28 of the First and Second Kinds</a> 29</h3></div></div></div> 30<h5> 31<a name="math_toolkit.bessel.mbessel.h0"></a> 32 <span class="phrase"><a name="math_toolkit.bessel.mbessel.synopsis"></a></span><a class="link" href="mbessel.html#math_toolkit.bessel.mbessel.synopsis">Synopsis</a> 33 </h5> 34<p> 35 <code class="computeroutput"><span class="preprocessor">#include</span> <span class="special"><</span><span class="identifier">boost</span><span class="special">/</span><span class="identifier">math</span><span class="special">/</span><span class="identifier">special_functions</span><span class="special">/</span><span class="identifier">bessel</span><span class="special">.</span><span class="identifier">hpp</span><span class="special">></span></code> 36 </p> 37<pre class="programlisting"><span class="keyword">template</span> <span class="special"><</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> 38<a class="link" href="../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_i</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> 39 40<span class="keyword">template</span> <span class="special"><</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="Chapter 21. Policies: Controlling Precision, Error Handling etc">Policy</a><span class="special">></span> 41<a class="link" href="../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_i</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="Chapter 21. Policies: Controlling Precision, Error Handling etc">Policy</a><span class="special">&);</span> 42 43<span class="keyword">template</span> <span class="special"><</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> 44<a class="link" href="../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_k</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> 45 46<span class="keyword">template</span> <span class="special"><</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="Chapter 21. Policies: Controlling Precision, Error Handling etc">Policy</a><span class="special">></span> 47<a class="link" href="../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_k</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="Chapter 21. Policies: Controlling Precision, Error Handling etc">Policy</a><span class="special">&);</span> 48</pre> 49<h5> 50<a name="math_toolkit.bessel.mbessel.h1"></a> 51 <span class="phrase"><a name="math_toolkit.bessel.mbessel.description"></a></span><a class="link" href="mbessel.html#math_toolkit.bessel.mbessel.description">Description</a> 52 </h5> 53<p> 54 The functions <a class="link" href="mbessel.html" title="Modified Bessel Functions of the First and Second Kinds">cyl_bessel_i</a> 55 and <a class="link" href="mbessel.html" title="Modified Bessel Functions of the First and Second Kinds">cyl_bessel_k</a> return 56 the result of the modified Bessel functions of the first and second kind 57 respectively: 58 </p> 59<div class="blockquote"><blockquote class="blockquote"><p> 60 cyl_bessel_i(v, x) = I<sub>v</sub>(x) 61 </p></blockquote></div> 62<div class="blockquote"><blockquote class="blockquote"><p> 63 cyl_bessel_k(v, x) = K<sub>v</sub>(x) 64 </p></blockquote></div> 65<p> 66 where: 67 </p> 68<div class="blockquote"><blockquote class="blockquote"><p> 69 <span class="inlinemediaobject"><img src="../../../equations/mbessel2.svg"></span> 70 71 </p></blockquote></div> 72<div class="blockquote"><blockquote class="blockquote"><p> 73 <span class="inlinemediaobject"><img src="../../../equations/mbessel3.svg"></span> 74 75 </p></blockquote></div> 76<p> 77 The return type of these functions is computed using the <a class="link" href="../result_type.html" title="Calculation of the Type of the Result"><span class="emphasis"><em>result 78 type calculation rules</em></span></a> when T1 and T2 are different types. 79 The functions are also optimised for the relatively common case that T1 is 80 an integer. 81 </p> 82<p> 83 The final <a class="link" href="../../policy.html" title="Chapter 21. Policies: Controlling Precision, Error Handling etc">Policy</a> argument is optional and can 84 be used to control the behaviour of the function: how it handles errors, 85 what level of precision to use etc. Refer to the <a class="link" href="../../policy.html" title="Chapter 21. Policies: Controlling Precision, Error Handling etc">policy 86 documentation for more details</a>. 87 </p> 88<p> 89 The functions return the result of <a class="link" href="../error_handling.html#math_toolkit.error_handling.domain_error">domain_error</a> 90 whenever the result is undefined or complex. For <a class="link" href="bessel_first.html" title="Bessel Functions of the First and Second Kinds">cyl_bessel_j</a> 91 this occurs when <code class="computeroutput"><span class="identifier">x</span> <span class="special"><</span> 92 <span class="number">0</span></code> and v is not an integer, or when 93 <code class="computeroutput"><span class="identifier">x</span> <span class="special">==</span> 94 <span class="number">0</span></code> and <code class="computeroutput"><span class="identifier">v</span> 95 <span class="special">!=</span> <span class="number">0</span></code>. 96 For <a class="link" href="bessel_first.html" title="Bessel Functions of the First and Second Kinds">cyl_neumann</a> this 97 occurs when <code class="computeroutput"><span class="identifier">x</span> <span class="special"><=</span> 98 <span class="number">0</span></code>. 99 </p> 100<p> 101 The following graph illustrates the exponential behaviour of I<sub>v</sub>. 102 </p> 103<div class="blockquote"><blockquote class="blockquote"><p> 104 <span class="inlinemediaobject"><img src="../../../graphs/cyl_bessel_i.svg" align="middle"></span> 105 106 </p></blockquote></div> 107<p> 108 The following graph illustrates the exponential decay of K<sub>v</sub>. 109 </p> 110<div class="blockquote"><blockquote class="blockquote"><p> 111 <span class="inlinemediaobject"><img src="../../../graphs/cyl_bessel_k.svg" align="middle"></span> 112 113 </p></blockquote></div> 114<h5> 115<a name="math_toolkit.bessel.mbessel.h2"></a> 116 <span class="phrase"><a name="math_toolkit.bessel.mbessel.testing"></a></span><a class="link" href="mbessel.html#math_toolkit.bessel.mbessel.testing">Testing</a> 117 </h5> 118<p> 119 There are two sets of test values: spot values calculated using <a href="http://functions.wolfram.com" target="_top">functions.wolfram.com</a>, 120 and a much larger set of tests computed using a simplified version of this 121 implementation (with all the special case handling removed). 122 </p> 123<h5> 124<a name="math_toolkit.bessel.mbessel.h3"></a> 125 <span class="phrase"><a name="math_toolkit.bessel.mbessel.accuracy"></a></span><a class="link" href="mbessel.html#math_toolkit.bessel.mbessel.accuracy">Accuracy</a> 126 </h5> 127<p> 128 The following tables show how the accuracy of these functions varies on various 129 platforms, along with comparison to other libraries. Note that only results 130 for the widest floating-point type on the system are given, as narrower types 131 have <a class="link" href="../relative_error.html#math_toolkit.relative_error.zero_error">effectively zero 132 error</a>. All values are relative errors in units of epsilon. Note that 133 our test suite includes some fairly extreme inputs which results in most 134 of the worst problem cases in other libraries: 135 </p> 136<div class="table"> 137<a name="math_toolkit.bessel.mbessel.table_cyl_bessel_i_integer_orders_"></a><p class="title"><b>Table 8.44. Error rates for cyl_bessel_i (integer orders)</b></p> 138<div class="table-contents"><table class="table" summary="Error rates for cyl_bessel_i (integer orders)"> 139<colgroup> 140<col> 141<col> 142<col> 143<col> 144<col> 145</colgroup> 146<thead><tr> 147<th> 148 </th> 149<th> 150 <p> 151 GNU C++ version 7.1.0<br> linux<br> double 152 </p> 153 </th> 154<th> 155 <p> 156 GNU C++ version 7.1.0<br> linux<br> long double 157 </p> 158 </th> 159<th> 160 <p> 161 Sun compiler version 0x5150<br> Sun Solaris<br> long double 162 </p> 163 </th> 164<th> 165 <p> 166 Microsoft Visual C++ version 14.1<br> Win32<br> double 167 </p> 168 </th> 169</tr></thead> 170<tbody> 171<tr> 172<td> 173 <p> 174 Bessel I0: Mathworld Data (Integer Version) 175 </p> 176 </td> 177<td> 178 <p> 179 <span class="blue">Max = 0ε (Mean = 0ε)</span><br> <br> (<span class="emphasis"><em>GSL 180 2.1:</em></span> Max = 0.79ε (Mean = 0.482ε))<br> (<span class="emphasis"><em>Rmath 181 3.2.3:</em></span> Max = 1.52ε (Mean = 0.622ε) <a class="link" href="../logs_and_tables/logs.html#errors_GNU_C_version_7_1_0_linux_double_cyl_bessel_i_integer_orders__Rmath_3_2_3_Bessel_I0_Mathworld_Data_Integer_Version_">And 182 other failures.</a>) 183 </p> 184 </td> 185<td> 186 <p> 187 <span class="blue">Max = 1.95ε (Mean = 0.738ε)</span><br> <br> 188 (<span class="emphasis"><em><cmath>:</em></span> Max = 8.49ε (Mean = 3.46ε) <a class="link" href="../logs_and_tables/logs.html#errors_GNU_C_version_7_1_0_linux_long_double_cyl_bessel_i_integer_orders___cmath__Bessel_I0_Mathworld_Data_Integer_Version_">And 189 other failures.</a>) 190 </p> 191 </td> 192<td> 193 <p> 194 <span class="blue">Max = 1.95ε (Mean = 0.661ε)</span> 195 </p> 196 </td> 197<td> 198 <p> 199 <span class="blue">Max = 0.762ε (Mean = 0.329ε)</span> 200 </p> 201 </td> 202</tr> 203<tr> 204<td> 205 <p> 206 Bessel I1: Mathworld Data (Integer Version) 207 </p> 208 </td> 209<td> 210 <p> 211 <span class="blue">Max = 0ε (Mean = 0ε)</span><br> <br> (<span class="emphasis"><em>GSL 212 2.1:</em></span> Max = 0.82ε (Mean = 0.456ε))<br> (<span class="emphasis"><em>Rmath 213 3.2.3:</em></span> Max = 1.53ε (Mean = 0.483ε) <a class="link" href="../logs_and_tables/logs.html#errors_GNU_C_version_7_1_0_linux_double_cyl_bessel_i_integer_orders__Rmath_3_2_3_Bessel_I1_Mathworld_Data_Integer_Version_">And 214 other failures.</a>) 215 </p> 216 </td> 217<td> 218 <p> 219 <span class="blue">Max = 0.64ε (Mean = 0.202ε)</span><br> <br> 220 (<span class="emphasis"><em><cmath>:</em></span> Max = 5ε (Mean = 2.15ε) <a class="link" href="../logs_and_tables/logs.html#errors_GNU_C_version_7_1_0_linux_long_double_cyl_bessel_i_integer_orders___cmath__Bessel_I1_Mathworld_Data_Integer_Version_">And 221 other failures.</a>) 222 </p> 223 </td> 224<td> 225 <p> 226 <span class="blue">Max = 0.64ε (Mean = 0.202ε)</span> 227 </p> 228 </td> 229<td> 230 <p> 231 <span class="blue">Max = 0.767ε (Mean = 0.398ε)</span> 232 </p> 233 </td> 234</tr> 235<tr> 236<td> 237 <p> 238 Bessel In: Mathworld Data (Integer Version) 239 </p> 240 </td> 241<td> 242 <p> 243 <span class="blue">Max = 0ε (Mean = 0ε)</span><br> <br> (<span class="emphasis"><em>GSL 244 2.1:</em></span> Max = 5.15ε (Mean = 2.13ε) <a class="link" href="../logs_and_tables/logs.html#errors_GNU_C_version_7_1_0_linux_double_cyl_bessel_i_integer_orders__GSL_2_1_Bessel_In_Mathworld_Data_Integer_Version_">And 245 other failures.</a>)<br> (<span class="emphasis"><em>Rmath 3.2.3:</em></span> 246 Max = 1.73ε (Mean = 0.601ε) <a class="link" href="../logs_and_tables/logs.html#errors_GNU_C_version_7_1_0_linux_double_cyl_bessel_i_integer_orders__Rmath_3_2_3_Bessel_In_Mathworld_Data_Integer_Version_">And 247 other failures.</a>) 248 </p> 249 </td> 250<td> 251 <p> 252 <span class="blue">Max = 1.8ε (Mean = 1.33ε)</span><br> <br> 253 (<span class="emphasis"><em><cmath>:</em></span> Max = 430ε (Mean = 163ε) <a class="link" href="../logs_and_tables/logs.html#errors_GNU_C_version_7_1_0_linux_long_double_cyl_bessel_i_integer_orders___cmath__Bessel_In_Mathworld_Data_Integer_Version_">And 254 other failures.</a>) 255 </p> 256 </td> 257<td> 258 <p> 259 <span class="blue">Max = 463ε (Mean = 140ε)</span> 260 </p> 261 </td> 262<td> 263 <p> 264 <span class="blue">Max = 3.46ε (Mean = 1.32ε)</span> 265 </p> 266 </td> 267</tr> 268</tbody> 269</table></div> 270</div> 271<br class="table-break"><div class="table"> 272<a name="math_toolkit.bessel.mbessel.table_cyl_bessel_i"></a><p class="title"><b>Table 8.45. Error rates for cyl_bessel_i</b></p> 273<div class="table-contents"><table class="table" summary="Error rates for cyl_bessel_i"> 274<colgroup> 275<col> 276<col> 277<col> 278<col> 279<col> 280</colgroup> 281<thead><tr> 282<th> 283 </th> 284<th> 285 <p> 286 GNU C++ version 7.1.0<br> linux<br> double 287 </p> 288 </th> 289<th> 290 <p> 291 GNU C++ version 7.1.0<br> linux<br> long double 292 </p> 293 </th> 294<th> 295 <p> 296 Sun compiler version 0x5150<br> Sun Solaris<br> long double 297 </p> 298 </th> 299<th> 300 <p> 301 Microsoft Visual C++ version 14.1<br> Win32<br> double 302 </p> 303 </th> 304</tr></thead> 305<tbody> 306<tr> 307<td> 308 <p> 309 Bessel I0: Mathworld Data 310 </p> 311 </td> 312<td> 313 <p> 314 <span class="blue">Max = 0ε (Mean = 0ε)</span><br> <br> (<span class="emphasis"><em>GSL 315 2.1:</em></span> Max = 270ε (Mean = 91.6ε) <a class="link" href="../logs_and_tables/logs.html#errors_GNU_C_version_7_1_0_linux_double_cyl_bessel_i_GSL_2_1_Bessel_I0_Mathworld_Data">And 316 other failures.</a>)<br> (<span class="emphasis"><em>Rmath 3.2.3:</em></span> 317 Max = 1.52ε (Mean = 0.622ε) <a class="link" href="../logs_and_tables/logs.html#errors_GNU_C_version_7_1_0_linux_double_cyl_bessel_i_Rmath_3_2_3_Bessel_I0_Mathworld_Data">And 318 other failures.</a>) 319 </p> 320 </td> 321<td> 322 <p> 323 <span class="blue">Max = 1.95ε (Mean = 0.738ε)</span><br> <br> 324 (<span class="emphasis"><em><cmath>:</em></span> Max = 8.49ε (Mean = 3.46ε) <a class="link" href="../logs_and_tables/logs.html#errors_GNU_C_version_7_1_0_linux_long_double_cyl_bessel_i__cmath__Bessel_I0_Mathworld_Data">And 325 other failures.</a>) 326 </p> 327 </td> 328<td> 329 <p> 330 <span class="blue">Max = 1.95ε (Mean = 0.661ε)</span> 331 </p> 332 </td> 333<td> 334 <p> 335 <span class="blue">Max = 0.762ε (Mean = 0.329ε)</span> 336 </p> 337 </td> 338</tr> 339<tr> 340<td> 341 <p> 342 Bessel I1: Mathworld Data 343 </p> 344 </td> 345<td> 346 <p> 347 <span class="blue">Max = 0ε (Mean = 0ε)</span><br> <br> (<span class="emphasis"><em>GSL 348 2.1:</em></span> Max = 128ε (Mean = 41ε) <a class="link" href="../logs_and_tables/logs.html#errors_GNU_C_version_7_1_0_linux_double_cyl_bessel_i_GSL_2_1_Bessel_I1_Mathworld_Data">And 349 other failures.</a>)<br> (<span class="emphasis"><em>Rmath 3.2.3:</em></span> 350 Max = 1.53ε (Mean = 0.483ε) <a class="link" href="../logs_and_tables/logs.html#errors_GNU_C_version_7_1_0_linux_double_cyl_bessel_i_Rmath_3_2_3_Bessel_I1_Mathworld_Data">And 351 other failures.</a>) 352 </p> 353 </td> 354<td> 355 <p> 356 <span class="blue">Max = 0.64ε (Mean = 0.202ε)</span><br> <br> 357 (<span class="emphasis"><em><cmath>:</em></span> Max = 5ε (Mean = 2.15ε) <a class="link" href="../logs_and_tables/logs.html#errors_GNU_C_version_7_1_0_linux_long_double_cyl_bessel_i__cmath__Bessel_I1_Mathworld_Data">And 358 other failures.</a>) 359 </p> 360 </td> 361<td> 362 <p> 363 <span class="blue">Max = 0.64ε (Mean = 0.202ε)</span> 364 </p> 365 </td> 366<td> 367 <p> 368 <span class="blue">Max = 0.767ε (Mean = 0.398ε)</span> 369 </p> 370 </td> 371</tr> 372<tr> 373<td> 374 <p> 375 Bessel In: Mathworld Data 376 </p> 377 </td> 378<td> 379 <p> 380 <span class="blue">Max = 0ε (Mean = 0ε)</span><br> <br> (<span class="emphasis"><em>GSL 381 2.1:</em></span> Max = 2.31ε (Mean = 0.838ε) <a class="link" href="../logs_and_tables/logs.html#errors_GNU_C_version_7_1_0_linux_double_cyl_bessel_i_GSL_2_1_Bessel_In_Mathworld_Data">And 382 other failures.</a>)<br> (<span class="emphasis"><em>Rmath 3.2.3:</em></span> 383 Max = 1.73ε (Mean = 0.601ε) <a class="link" href="../logs_and_tables/logs.html#errors_GNU_C_version_7_1_0_linux_double_cyl_bessel_i_Rmath_3_2_3_Bessel_In_Mathworld_Data">And 384 other failures.</a>) 385 </p> 386 </td> 387<td> 388 <p> 389 <span class="blue">Max = 1.8ε (Mean = 1.33ε)</span><br> <br> 390 (<span class="emphasis"><em><cmath>:</em></span> Max = 430ε (Mean = 163ε) <a class="link" href="../logs_and_tables/logs.html#errors_GNU_C_version_7_1_0_linux_long_double_cyl_bessel_i__cmath__Bessel_In_Mathworld_Data">And 391 other failures.</a>) 392 </p> 393 </td> 394<td> 395 <p> 396 <span class="blue">Max = 463ε (Mean = 140ε)</span> 397 </p> 398 </td> 399<td> 400 <p> 401 <span class="blue">Max = 3.46ε (Mean = 1.32ε)</span> 402 </p> 403 </td> 404</tr> 405<tr> 406<td> 407 <p> 408 Bessel Iv: Mathworld Data 409 </p> 410 </td> 411<td> 412 <p> 413 <span class="blue">Max = 0ε (Mean = 0ε)</span><br> <br> (<span class="emphasis"><em>GSL 414 2.1:</em></span> Max = 5.95ε (Mean = 2.08ε) <a class="link" href="../logs_and_tables/logs.html#errors_GNU_C_version_7_1_0_linux_double_cyl_bessel_i_GSL_2_1_Bessel_Iv_Mathworld_Data">And 415 other failures.</a>)<br> (<span class="emphasis"><em>Rmath 3.2.3:</em></span> 416 Max = 3.53ε (Mean = 1.39ε)) 417 </p> 418 </td> 419<td> 420 <p> 421 <span class="blue">Max = 4.12ε (Mean = 1.85ε)</span><br> <br> 422 (<span class="emphasis"><em><cmath>:</em></span> Max = 616ε (Mean = 221ε) <a class="link" href="../logs_and_tables/logs.html#errors_GNU_C_version_7_1_0_linux_long_double_cyl_bessel_i__cmath__Bessel_Iv_Mathworld_Data">And 423 other failures.</a>) 424 </p> 425 </td> 426<td> 427 <p> 428 <span class="blue">Max = 4.12ε (Mean = 1.95ε)</span> 429 </p> 430 </td> 431<td> 432 <p> 433 <span class="blue">Max = 2.97ε (Mean = 1.24ε)</span> 434 </p> 435 </td> 436</tr> 437<tr> 438<td> 439 <p> 440 Bessel In: Random Data 441 </p> 442 </td> 443<td> 444 <p> 445 <span class="blue">Max = 0ε (Mean = 0ε)</span><br> <br> (<span class="emphasis"><em>GSL 446 2.1:</em></span> Max = 261ε (Mean = 53.2ε) <a class="link" href="../logs_and_tables/logs.html#errors_GNU_C_version_7_1_0_linux_double_cyl_bessel_i_GSL_2_1_Bessel_In_Random_Data">And 447 other failures.</a>)<br> (<span class="emphasis"><em>Rmath 3.2.3:</em></span> 448 Max = 7.37ε (Mean = 2.4ε)) 449 </p> 450 </td> 451<td> 452 <p> 453 <span class="blue">Max = 4.62ε (Mean = 1.06ε)</span><br> <br> 454 (<span class="emphasis"><em><cmath>:</em></span> Max = 645ε (Mean = 132ε)) 455 </p> 456 </td> 457<td> 458 <p> 459 <span class="blue">Max = 176ε (Mean = 39.1ε)</span> 460 </p> 461 </td> 462<td> 463 <p> 464 <span class="blue">Max = 9.67ε (Mean = 1.88ε)</span> 465 </p> 466 </td> 467</tr> 468<tr> 469<td> 470 <p> 471 Bessel Iv: Random Data 472 </p> 473 </td> 474<td> 475 <p> 476 <span class="blue">Max = 0.661ε (Mean = 0.0441ε)</span><br> 477 <br> (<span class="emphasis"><em>GSL 2.1:</em></span> Max = 6.18e+03ε (Mean = 1.55e+03ε) 478 <a class="link" href="../logs_and_tables/logs.html#errors_GNU_C_version_7_1_0_linux_double_cyl_bessel_i_GSL_2_1_Bessel_Iv_Random_Data">And 479 other failures.</a>)<br> (<span class="emphasis"><em>Rmath 3.2.3:</em></span> 480 <span class="red">Max = 4.28e+08ε (Mean = 2.85e+07ε))</span> 481 </p> 482 </td> 483<td> 484 <p> 485 <span class="blue">Max = 8.35ε (Mean = 1.62ε)</span><br> <br> 486 (<span class="emphasis"><em><cmath>:</em></span> Max = 1.05e+03ε (Mean = 224ε) 487 <a class="link" href="../logs_and_tables/logs.html#errors_GNU_C_version_7_1_0_linux_long_double_cyl_bessel_i__cmath__Bessel_Iv_Random_Data">And 488 other failures.</a>) 489 </p> 490 </td> 491<td> 492 <p> 493 <span class="blue">Max = 283ε (Mean = 88.4ε)</span> 494 </p> 495 </td> 496<td> 497 <p> 498 <span class="blue">Max = 7.46ε (Mean = 1.71ε)</span> 499 </p> 500 </td> 501</tr> 502<tr> 503<td> 504 <p> 505 Bessel Iv: Mathworld Data (large values) 506 </p> 507 </td> 508<td> 509 <p> 510 <span class="blue">Max = 0ε (Mean = 0ε)</span><br> <br> (<span class="emphasis"><em>GSL 511 2.1:</em></span> Max = 37ε (Mean = 18ε) <a class="link" href="../logs_and_tables/logs.html#errors_GNU_C_version_7_1_0_linux_double_cyl_bessel_i_GSL_2_1_Bessel_Iv_Mathworld_Data_large_values_">And 512 other failures.</a>)<br> (<span class="emphasis"><em>Rmath 3.2.3:</em></span> 513 <span class="red">Max = 3.77e+168ε (Mean = 2.39e+168ε) <a class="link" href="../logs_and_tables/logs.html#errors_GNU_C_version_7_1_0_linux_double_cyl_bessel_i_Rmath_3_2_3_Bessel_Iv_Mathworld_Data_large_values_">And 514 other failures.</a>)</span> 515 </p> 516 </td> 517<td> 518 <p> 519 <span class="blue">Max = 14.7ε (Mean = 6.66ε)</span><br> <br> 520 (<span class="emphasis"><em><cmath>:</em></span> Max = 118ε (Mean = 57.2ε) <a class="link" href="../logs_and_tables/logs.html#errors_GNU_C_version_7_1_0_linux_long_double_cyl_bessel_i__cmath__Bessel_Iv_Mathworld_Data_large_values_">And 521 other failures.</a>) 522 </p> 523 </td> 524<td> 525 <p> 526 <span class="blue">Max = 14.7ε (Mean = 6.59ε)</span> 527 </p> 528 </td> 529<td> 530 <p> 531 <span class="blue">Max = 3.67ε (Mean = 1.64ε)</span> 532 </p> 533 </td> 534</tr> 535</tbody> 536</table></div> 537</div> 538<br class="table-break"><div class="table"> 539<a name="math_toolkit.bessel.mbessel.table_cyl_bessel_k_integer_orders_"></a><p class="title"><b>Table 8.46. Error rates for cyl_bessel_k (integer orders)</b></p> 540<div class="table-contents"><table class="table" summary="Error rates for cyl_bessel_k (integer orders)"> 541<colgroup> 542<col> 543<col> 544<col> 545<col> 546<col> 547</colgroup> 548<thead><tr> 549<th> 550 </th> 551<th> 552 <p> 553 GNU C++ version 7.1.0<br> linux<br> long double 554 </p> 555 </th> 556<th> 557 <p> 558 GNU C++ version 7.1.0<br> linux<br> double 559 </p> 560 </th> 561<th> 562 <p> 563 Sun compiler version 0x5150<br> Sun Solaris<br> long double 564 </p> 565 </th> 566<th> 567 <p> 568 Microsoft Visual C++ version 14.1<br> Win32<br> double 569 </p> 570 </th> 571</tr></thead> 572<tbody> 573<tr> 574<td> 575 <p> 576 Bessel K0: Mathworld Data (Integer Version) 577 </p> 578 </td> 579<td> 580 <p> 581 <span class="blue">Max = 0.833ε (Mean = 0.436ε)</span><br> <br> 582 (<span class="emphasis"><em><cmath>:</em></span> Max = 9.33ε (Mean = 3.25ε)) 583 </p> 584 </td> 585<td> 586 <p> 587 <span class="blue">Max = 0ε (Mean = 0ε)</span><br> <br> (<span class="emphasis"><em>GSL 588 2.1:</em></span> Max = 1.2ε (Mean = 0.733ε))<br> (<span class="emphasis"><em>Rmath 589 3.2.3:</em></span> Max = 0.833ε (Mean = 0.601ε)) 590 </p> 591 </td> 592<td> 593 <p> 594 <span class="blue">Max = 0.833ε (Mean = 0.436ε)</span> 595 </p> 596 </td> 597<td> 598 <p> 599 <span class="blue">Max = 0.833ε (Mean = 0.552ε)</span> 600 </p> 601 </td> 602</tr> 603<tr> 604<td> 605 <p> 606 Bessel K1: Mathworld Data (Integer Version) 607 </p> 608 </td> 609<td> 610 <p> 611 <span class="blue">Max = 0.786ε (Mean = 0.329ε)</span><br> <br> 612 (<span class="emphasis"><em><cmath>:</em></span> Max = 8.94ε (Mean = 3.19ε)) 613 </p> 614 </td> 615<td> 616 <p> 617 <span class="blue">Max = 0ε (Mean = 0ε)</span><br> <br> (<span class="emphasis"><em>GSL 618 2.1:</em></span> Max = 0.626ε (Mean = 0.333ε))<br> (<span class="emphasis"><em>Rmath 619 3.2.3:</em></span> Max = 0.894ε (Mean = 0.516ε)) 620 </p> 621 </td> 622<td> 623 <p> 624 <span class="blue">Max = 0.786ε (Mean = 0.329ε)</span> 625 </p> 626 </td> 627<td> 628 <p> 629 <span class="blue">Max = 0.786ε (Mean = 0.39ε)</span> 630 </p> 631 </td> 632</tr> 633<tr> 634<td> 635 <p> 636 Bessel Kn: Mathworld Data (Integer Version) 637 </p> 638 </td> 639<td> 640 <p> 641 <span class="blue">Max = 2.6ε (Mean = 1.21ε)</span><br> <br> 642 (<span class="emphasis"><em><cmath>:</em></span> Max = 12.9ε (Mean = 4.91ε) <a class="link" href="../logs_and_tables/logs.html#errors_GNU_C_version_7_1_0_linux_long_double_cyl_bessel_k_integer_orders___cmath__Bessel_Kn_Mathworld_Data_Integer_Version_">And 643 other failures.</a>) 644 </p> 645 </td> 646<td> 647 <p> 648 <span class="blue">Max = 0ε (Mean = 0ε)</span><br> <br> (<span class="emphasis"><em>GSL 649 2.1:</em></span> Max = 168ε (Mean = 59.5ε))<br> (<span class="emphasis"><em>Rmath 650 3.2.3:</em></span> Max = 8.48ε (Mean = 2.98ε)) 651 </p> 652 </td> 653<td> 654 <p> 655 <span class="blue">Max = 2.6ε (Mean = 1.21ε)</span> 656 </p> 657 </td> 658<td> 659 <p> 660 <span class="blue">Max = 3.63ε (Mean = 1.46ε)</span> 661 </p> 662 </td> 663</tr> 664</tbody> 665</table></div> 666</div> 667<br class="table-break"><div class="table"> 668<a name="math_toolkit.bessel.mbessel.table_cyl_bessel_k"></a><p class="title"><b>Table 8.47. Error rates for cyl_bessel_k</b></p> 669<div class="table-contents"><table class="table" summary="Error rates for cyl_bessel_k"> 670<colgroup> 671<col> 672<col> 673<col> 674<col> 675<col> 676</colgroup> 677<thead><tr> 678<th> 679 </th> 680<th> 681 <p> 682 GNU C++ version 7.1.0<br> linux<br> long double 683 </p> 684 </th> 685<th> 686 <p> 687 GNU C++ version 7.1.0<br> linux<br> double 688 </p> 689 </th> 690<th> 691 <p> 692 Sun compiler version 0x5150<br> Sun Solaris<br> long double 693 </p> 694 </th> 695<th> 696 <p> 697 Microsoft Visual C++ version 14.1<br> Win32<br> double 698 </p> 699 </th> 700</tr></thead> 701<tbody> 702<tr> 703<td> 704 <p> 705 Bessel K0: Mathworld Data 706 </p> 707 </td> 708<td> 709 <p> 710 <span class="blue">Max = 0.833ε (Mean = 0.436ε)</span><br> <br> 711 (<span class="emphasis"><em><cmath>:</em></span> Max = 9.33ε (Mean = 3.25ε)) 712 </p> 713 </td> 714<td> 715 <p> 716 <span class="blue">Max = 0ε (Mean = 0ε)</span><br> <br> (<span class="emphasis"><em>GSL 717 2.1:</em></span> Max = 6.04ε (Mean = 2.16ε))<br> (<span class="emphasis"><em>Rmath 718 3.2.3:</em></span> Max = 0.833ε (Mean = 0.601ε)) 719 </p> 720 </td> 721<td> 722 <p> 723 <span class="blue">Max = 0.833ε (Mean = 0.436ε)</span> 724 </p> 725 </td> 726<td> 727 <p> 728 <span class="blue">Max = 0.833ε (Mean = 0.552ε)</span> 729 </p> 730 </td> 731</tr> 732<tr> 733<td> 734 <p> 735 Bessel K1: Mathworld Data 736 </p> 737 </td> 738<td> 739 <p> 740 <span class="blue">Max = 0.786ε (Mean = 0.329ε)</span><br> <br> 741 (<span class="emphasis"><em><cmath>:</em></span> Max = 8.94ε (Mean = 3.19ε)) 742 </p> 743 </td> 744<td> 745 <p> 746 <span class="blue">Max = 0ε (Mean = 0ε)</span><br> <br> (<span class="emphasis"><em>GSL 747 2.1:</em></span> Max = 6.26ε (Mean = 2.21ε))<br> (<span class="emphasis"><em>Rmath 748 3.2.3:</em></span> Max = 0.894ε (Mean = 0.516ε)) 749 </p> 750 </td> 751<td> 752 <p> 753 <span class="blue">Max = 0.786ε (Mean = 0.329ε)</span> 754 </p> 755 </td> 756<td> 757 <p> 758 <span class="blue">Max = 0.786ε (Mean = 0.39ε)</span> 759 </p> 760 </td> 761</tr> 762<tr> 763<td> 764 <p> 765 Bessel Kn: Mathworld Data 766 </p> 767 </td> 768<td> 769 <p> 770 <span class="blue">Max = 2.6ε (Mean = 1.21ε)</span><br> <br> 771 (<span class="emphasis"><em><cmath>:</em></span> Max = 12.9ε (Mean = 4.91ε) <a class="link" href="../logs_and_tables/logs.html#errors_GNU_C_version_7_1_0_linux_long_double_cyl_bessel_k__cmath__Bessel_Kn_Mathworld_Data">And 772 other failures.</a>) 773 </p> 774 </td> 775<td> 776 <p> 777 <span class="blue">Max = 0ε (Mean = 0ε)</span><br> <br> (<span class="emphasis"><em>GSL 778 2.1:</em></span> Max = 3.36ε (Mean = 1.43ε) <a class="link" href="../logs_and_tables/logs.html#errors_GNU_C_version_7_1_0_linux_double_cyl_bessel_k_GSL_2_1_Bessel_Kn_Mathworld_Data">And 779 other failures.</a>)<br> (<span class="emphasis"><em>Rmath 3.2.3:</em></span> 780 Max = 8.48ε (Mean = 2.98ε)) 781 </p> 782 </td> 783<td> 784 <p> 785 <span class="blue">Max = 2.6ε (Mean = 1.21ε)</span> 786 </p> 787 </td> 788<td> 789 <p> 790 <span class="blue">Max = 3.63ε (Mean = 1.46ε)</span> 791 </p> 792 </td> 793</tr> 794<tr> 795<td> 796 <p> 797 Bessel Kv: Mathworld Data 798 </p> 799 </td> 800<td> 801 <p> 802 <span class="blue">Max = 3.58ε (Mean = 2.39ε)</span><br> <br> 803 (<span class="emphasis"><em><cmath>:</em></span> Max = 13ε (Mean = 4.81ε) <a class="link" href="../logs_and_tables/logs.html#errors_GNU_C_version_7_1_0_linux_long_double_cyl_bessel_k__cmath__Bessel_Kv_Mathworld_Data">And 804 other failures.</a>) 805 </p> 806 </td> 807<td> 808 <p> 809 <span class="blue">Max = 0ε (Mean = 0ε)</span><br> <br> (<span class="emphasis"><em>GSL 810 2.1:</em></span> Max = 5.47ε (Mean = 2.04ε) <a class="link" href="../logs_and_tables/logs.html#errors_GNU_C_version_7_1_0_linux_double_cyl_bessel_k_GSL_2_1_Bessel_Kv_Mathworld_Data">And 811 other failures.</a>)<br> (<span class="emphasis"><em>Rmath 3.2.3:</em></span> 812 Max = 3.15ε (Mean = 1.35ε)) 813 </p> 814 </td> 815<td> 816 <p> 817 <span class="blue">Max = 5.21ε (Mean = 2.53ε)</span> 818 </p> 819 </td> 820<td> 821 <p> 822 <span class="blue">Max = 4.78ε (Mean = 2.19ε)</span> 823 </p> 824 </td> 825</tr> 826<tr> 827<td> 828 <p> 829 Bessel Kv: Mathworld Data (large values) 830 </p> 831 </td> 832<td> 833 <p> 834 <span class="blue">Max = 42.3ε (Mean = 21ε)</span><br> <br> 835 (<span class="emphasis"><em><cmath>:</em></span> Max = 42.3ε (Mean = 19.8ε) <a class="link" href="../logs_and_tables/logs.html#errors_GNU_C_version_7_1_0_linux_long_double_cyl_bessel_k__cmath__Bessel_Kv_Mathworld_Data_large_values_">And 836 other failures.</a>) 837 </p> 838 </td> 839<td> 840 <p> 841 <span class="blue">Max = 0ε (Mean = 0ε)</span><br> <br> (<span class="emphasis"><em>GSL 842 2.1:</em></span> Max = 308ε (Mean = 142ε) <a class="link" href="../logs_and_tables/logs.html#errors_GNU_C_version_7_1_0_linux_double_cyl_bessel_k_GSL_2_1_Bessel_Kv_Mathworld_Data_large_values_">And 843 other failures.</a>)<br> (<span class="emphasis"><em>Rmath 3.2.3:</em></span> 844 Max = 84.6ε (Mean = 37.8ε)) 845 </p> 846 </td> 847<td> 848 <p> 849 <span class="blue">Max = 42.3ε (Mean = 21ε)</span> 850 </p> 851 </td> 852<td> 853 <p> 854 <span class="blue">Max = 59.8ε (Mean = 26.9ε)</span> 855 </p> 856 </td> 857</tr> 858<tr> 859<td> 860 <p> 861 Bessel Kn: Random Data 862 </p> 863 </td> 864<td> 865 <p> 866 <span class="blue">Max = 4.55ε (Mean = 1.12ε)</span><br> <br> 867 (<span class="emphasis"><em><cmath>:</em></span> Max = 13.9ε (Mean = 2.91ε)) 868 </p> 869 </td> 870<td> 871 <p> 872 <span class="blue">Max = 0.764ε (Mean = 0.0348ε)</span><br> 873 <br> (<span class="emphasis"><em>GSL 2.1:</em></span> Max = 8.71ε (Mean = 1.76ε) 874 <a class="link" href="../logs_and_tables/logs.html#errors_GNU_C_version_7_1_0_linux_double_cyl_bessel_k_GSL_2_1_Bessel_Kn_Random_Data">And 875 other failures.</a>)<br> (<span class="emphasis"><em>Rmath 3.2.3:</em></span> 876 Max = 7.47ε (Mean = 1.34ε)) 877 </p> 878 </td> 879<td> 880 <p> 881 <span class="blue">Max = 4.55ε (Mean = 1.12ε)</span> 882 </p> 883 </td> 884<td> 885 <p> 886 <span class="blue">Max = 9.34ε (Mean = 1.7ε)</span> 887 </p> 888 </td> 889</tr> 890<tr> 891<td> 892 <p> 893 Bessel Kv: Random Data 894 </p> 895 </td> 896<td> 897 <p> 898 <span class="blue">Max = 7.88ε (Mean = 1.48ε)</span><br> <br> 899 (<span class="emphasis"><em><cmath>:</em></span> Max = 13.6ε (Mean = 2.68ε) <a class="link" href="../logs_and_tables/logs.html#errors_GNU_C_version_7_1_0_linux_long_double_cyl_bessel_k__cmath__Bessel_Kv_Random_Data">And 900 other failures.</a>) 901 </p> 902 </td> 903<td> 904 <p> 905 <span class="blue">Max = 0.507ε (Mean = 0.0313ε)</span><br> 906 <br> (<span class="emphasis"><em>GSL 2.1:</em></span> Max = 9.71ε (Mean = 1.47ε) 907 <a class="link" href="../logs_and_tables/logs.html#errors_GNU_C_version_7_1_0_linux_double_cyl_bessel_k_GSL_2_1_Bessel_Kv_Random_Data">And 908 other failures.</a>)<br> (<span class="emphasis"><em>Rmath 3.2.3:</em></span> 909 Max = 7.37ε (Mean = 1.49ε)) 910 </p> 911 </td> 912<td> 913 <p> 914 <span class="blue">Max = 7.88ε (Mean = 1.47ε)</span> 915 </p> 916 </td> 917<td> 918 <p> 919 <span class="blue">Max = 8.33ε (Mean = 1.62ε)</span> 920 </p> 921 </td> 922</tr> 923</tbody> 924</table></div> 925</div> 926<br class="table-break"><p> 927 The following error plot are based on an exhaustive search of the functions 928 domain for I0, I1, K0, and K1, MSVC-15.5 at <code class="computeroutput"><span class="keyword">double</span></code> 929 precision, and GCC-7.1/Ubuntu for <code class="computeroutput"><span class="keyword">long</span> 930 <span class="keyword">double</span></code> and <code class="computeroutput"><span class="identifier">__float128</span></code>. 931 </p> 932<div class="blockquote"><blockquote class="blockquote"><p> 933 <span class="inlinemediaobject"><img src="../../../graphs/i0__double.svg" align="middle"></span> 934 935 </p></blockquote></div> 936<div class="blockquote"><blockquote class="blockquote"><p> 937 <span class="inlinemediaobject"><img src="../../../graphs/i0__80_bit_long_double.svg" align="middle"></span> 938 939 </p></blockquote></div> 940<div class="blockquote"><blockquote class="blockquote"><p> 941 <span class="inlinemediaobject"><img src="../../../graphs/i0____float128.svg" align="middle"></span> 942 943 </p></blockquote></div> 944<div class="blockquote"><blockquote class="blockquote"><p> 945 <span class="inlinemediaobject"><img src="../../../graphs/i1__double.svg" align="middle"></span> 946 947 </p></blockquote></div> 948<div class="blockquote"><blockquote class="blockquote"><p> 949 <span class="inlinemediaobject"><img src="../../../graphs/i1__80_bit_long_double.svg" align="middle"></span> 950 951 </p></blockquote></div> 952<div class="blockquote"><blockquote class="blockquote"><p> 953 <span class="inlinemediaobject"><img src="../../../graphs/i1____float128.svg" align="middle"></span> 954 955 </p></blockquote></div> 956<div class="blockquote"><blockquote class="blockquote"><p> 957 <span class="inlinemediaobject"><img src="../../../graphs/k0__double.svg" align="middle"></span> 958 959 </p></blockquote></div> 960<div class="blockquote"><blockquote class="blockquote"><p> 961 <span class="inlinemediaobject"><img src="../../../graphs/k0__80_bit_long_double.svg" align="middle"></span> 962 963 </p></blockquote></div> 964<div class="blockquote"><blockquote class="blockquote"><p> 965 <span class="inlinemediaobject"><img src="../../../graphs/k0____float128.svg" align="middle"></span> 966 967 </p></blockquote></div> 968<div class="blockquote"><blockquote class="blockquote"><p> 969 <span class="inlinemediaobject"><img src="../../../graphs/k1__double.svg" align="middle"></span> 970 971 </p></blockquote></div> 972<div class="blockquote"><blockquote class="blockquote"><p> 973 <span class="inlinemediaobject"><img src="../../../graphs/k1__80_bit_long_double.svg" align="middle"></span> 974 975 </p></blockquote></div> 976<div class="blockquote"><blockquote class="blockquote"><p> 977 <span class="inlinemediaobject"><img src="../../../graphs/k1____float128.svg" align="middle"></span> 978 979 </p></blockquote></div> 980<h5> 981<a name="math_toolkit.bessel.mbessel.h4"></a> 982 <span class="phrase"><a name="math_toolkit.bessel.mbessel.implementation"></a></span><a class="link" href="mbessel.html#math_toolkit.bessel.mbessel.implementation">Implementation</a> 983 </h5> 984<p> 985 The following are handled as special cases first: 986 </p> 987<p> 988 When computing I<sub>v</sub> for <span class="emphasis"><em>x < 0</em></span>, then ν must be an integer 989 or a domain error occurs. If ν is an integer, then the function is odd if ν is 990 odd and even if ν is even, and we can reflect to <span class="emphasis"><em>x > 0</em></span>. 991 </p> 992<p> 993 For I<sub>v</sub> with v equal to 0, 1 or 0.5 are handled as special cases. 994 </p> 995<p> 996 The 0 and 1 cases use polynomial approximations on finite and infinite intervals. 997 The approximating forms are based on <a href="http://www.advanpix.com/2015/11/11/rational-approximations-for-the-modified-bessel-function-of-the-first-kind-i0-computations-double-precision/" target="_top">"Rational 998 Approximations for the Modified Bessel Function of the First Kind - I<sub>0</sub>(x) 999 for Computations with Double Precision"</a> by Pavel Holoborodko, 1000 extended by us to deal with up to 128-bit precision (with different approximations 1001 for each target precision). 1002 </p> 1003<div class="blockquote"><blockquote class="blockquote"><p> 1004 <span class="inlinemediaobject"><img src="../../../equations/bessel21.svg"></span> 1005 1006 </p></blockquote></div> 1007<div class="blockquote"><blockquote class="blockquote"><p> 1008 <span class="inlinemediaobject"><img src="../../../equations/bessel20.svg"></span> 1009 1010 </p></blockquote></div> 1011<div class="blockquote"><blockquote class="blockquote"><p> 1012 <span class="inlinemediaobject"><img src="../../../equations/bessel17.svg"></span> 1013 1014 </p></blockquote></div> 1015<div class="blockquote"><blockquote class="blockquote"><p> 1016 <span class="inlinemediaobject"><img src="../../../equations/bessel18.svg"></span> 1017 1018 </p></blockquote></div> 1019<p> 1020 Similarly we have: 1021 </p> 1022<div class="blockquote"><blockquote class="blockquote"><p> 1023 <span class="inlinemediaobject"><img src="../../../equations/bessel22.svg"></span> 1024 1025 </p></blockquote></div> 1026<div class="blockquote"><blockquote class="blockquote"><p> 1027 <span class="inlinemediaobject"><img src="../../../equations/bessel23.svg"></span> 1028 1029 </p></blockquote></div> 1030<div class="blockquote"><blockquote class="blockquote"><p> 1031 <span class="inlinemediaobject"><img src="../../../equations/bessel24.svg"></span> 1032 1033 </p></blockquote></div> 1034<div class="blockquote"><blockquote class="blockquote"><p> 1035 <span class="inlinemediaobject"><img src="../../../equations/bessel25.svg"></span> 1036 1037 </p></blockquote></div> 1038<p> 1039 The 0.5 case is a simple trigonometric function: 1040 </p> 1041<div class="blockquote"><blockquote class="blockquote"><p> 1042 I<sub>0.5</sub>(x) = sqrt(2 / πx) * sinh(x) 1043 </p></blockquote></div> 1044<p> 1045 For K<sub>v</sub> with <span class="emphasis"><em>v</em></span> an integer, the result is calculated using 1046 the recurrence relation: 1047 </p> 1048<div class="blockquote"><blockquote class="blockquote"><p> 1049 <span class="inlinemediaobject"><img src="../../../equations/mbessel5.svg"></span> 1050 1051 </p></blockquote></div> 1052<p> 1053 starting from K<sub>0</sub> and K<sub>1</sub> which are calculated using rational the approximations 1054 above. These rational approximations are accurate to around 19 digits, and 1055 are therefore only used when T has no more than 64 binary digits of precision. 1056 </p> 1057<p> 1058 When <span class="emphasis"><em>x</em></span> is small compared to <span class="emphasis"><em>v</em></span>, 1059 I<sub>v</sub>x is best computed directly from the series: 1060 </p> 1061<div class="blockquote"><blockquote class="blockquote"><p> 1062 <span class="inlinemediaobject"><img src="../../../equations/mbessel17.svg"></span> 1063 1064 </p></blockquote></div> 1065<p> 1066 In the general case, we first normalize ν to [<code class="literal">0, [inf]</code>) 1067 with the help of the reflection formulae: 1068 </p> 1069<div class="blockquote"><blockquote class="blockquote"><p> 1070 <span class="inlinemediaobject"><img src="../../../equations/mbessel9.svg"></span> 1071 1072 </p></blockquote></div> 1073<div class="blockquote"><blockquote class="blockquote"><p> 1074 <span class="inlinemediaobject"><img src="../../../equations/mbessel10.svg"></span> 1075 1076 </p></blockquote></div> 1077<p> 1078 Let μ = ν - floor(ν + 1/2), then μ is the fractional part of ν such that |μ| <= 1/2 1079 (we need this for convergence later). The idea is to calculate K<sub>μ</sub>(x) and K<sub>μ+1</sub>(x), 1080 and use them to obtain I<sub>ν</sub>(x) and K<sub>ν</sub>(x). 1081 </p> 1082<p> 1083 The algorithm is proposed by Temme in 1084 </p> 1085<div class="blockquote"><blockquote class="blockquote"><p> 1086 N.M. Temme, <span class="emphasis"><em>On the numerical evaluation of the modified bessel 1087 function of the third kind</em></span>, Journal of Computational Physics, 1088 vol 19, 324 (1975), 1089 </p></blockquote></div> 1090<p> 1091 which needs two continued fractions as well as the Wronskian: 1092 </p> 1093<div class="blockquote"><blockquote class="blockquote"><p> 1094 <span class="inlinemediaobject"><img src="../../../equations/mbessel11.svg"></span> 1095 1096 </p></blockquote></div> 1097<div class="blockquote"><blockquote class="blockquote"><p> 1098 <span class="inlinemediaobject"><img src="../../../equations/mbessel12.svg"></span> 1099 1100 </p></blockquote></div> 1101<div class="blockquote"><blockquote class="blockquote"><p> 1102 <span class="inlinemediaobject"><img src="../../../equations/mbessel8.svg"></span> 1103 1104 </p></blockquote></div> 1105<p> 1106 The continued fractions are computed using the modified Lentz's method 1107 </p> 1108<div class="blockquote"><blockquote class="blockquote"><p> 1109 (W.J. Lentz, <span class="emphasis"><em>Generating Bessel functions in Mie scattering calculations 1110 using continued fractions</em></span>, Applied Optics, vol 15, 668 (1976)). 1111 </p></blockquote></div> 1112<p> 1113 Their convergence rates depend on <span class="emphasis"><em>x</em></span>, therefore we need 1114 different strategies for large <span class="emphasis"><em>x</em></span> and small <span class="emphasis"><em>x</em></span>. 1115 </p> 1116<p> 1117 <span class="emphasis"><em>x > v</em></span>, CF1 needs O(<span class="emphasis"><em>x</em></span>) iterations 1118 to converge, CF2 converges rapidly. 1119 </p> 1120<p> 1121 <span class="emphasis"><em>x <= v</em></span>, CF1 converges rapidly, CF2 fails to converge 1122 when <span class="emphasis"><em>x</em></span> <code class="literal">-></code> 0. 1123 </p> 1124<p> 1125 When <span class="emphasis"><em>x</em></span> is large (<span class="emphasis"><em>x</em></span> > 2), both 1126 continued fractions converge (CF1 may be slow for really large <span class="emphasis"><em>x</em></span>). 1127 K<sub>μ</sub> and K<sub>μ+1</sub> 1128can be calculated by 1129 </p> 1130<div class="blockquote"><blockquote class="blockquote"><p> 1131 <span class="inlinemediaobject"><img src="../../../equations/mbessel13.svg"></span> 1132 1133 </p></blockquote></div> 1134<p> 1135 where 1136 </p> 1137<div class="blockquote"><blockquote class="blockquote"><p> 1138 <span class="inlinemediaobject"><img src="../../../equations/mbessel14.svg"></span> 1139 1140 </p></blockquote></div> 1141<p> 1142 <span class="emphasis"><em>S</em></span> is also a series that is summed along with CF2, see 1143 </p> 1144<div class="blockquote"><blockquote class="blockquote"><p> 1145 I.J. Thompson and A.R. Barnett, <span class="emphasis"><em>Modified Bessel functions I_v 1146 and K_v of real order and complex argument to selected accuracy</em></span>, 1147 Computer Physics Communications, vol 47, 245 (1987). 1148 </p></blockquote></div> 1149<p> 1150 When <span class="emphasis"><em>x</em></span> is small (<span class="emphasis"><em>x</em></span> <= 2), CF2 1151 convergence may fail (but CF1 works very well). The solution here is Temme's 1152 series: 1153 </p> 1154<div class="blockquote"><blockquote class="blockquote"><p> 1155 <span class="inlinemediaobject"><img src="../../../equations/mbessel15.svg"></span> 1156 1157 </p></blockquote></div> 1158<p> 1159 where 1160 </p> 1161<div class="blockquote"><blockquote class="blockquote"><p> 1162 <span class="inlinemediaobject"><img src="../../../equations/mbessel16.svg"></span> 1163 1164 </p></blockquote></div> 1165<p> 1166 f<sub>k</sub> and h<sub>k</sub> 1167are also computed by recursions (involving gamma functions), but 1168 the formulas are a little complicated, readers are referred to 1169 </p> 1170<div class="blockquote"><blockquote class="blockquote"><p> 1171 N.M. Temme, <span class="emphasis"><em>On the numerical evaluation of the modified Bessel 1172 function of the third kind</em></span>, Journal of Computational Physics, 1173 vol 19, 324 (1975). 1174 </p></blockquote></div> 1175<p> 1176 Note: Temme's series converge only for |μ| <= 1/2. 1177 </p> 1178<p> 1179 K<sub>ν</sub>(x) is then calculated from the forward recurrence, as is K<sub>ν+1</sub>(x). With these 1180 two values and f<sub>ν</sub>, the Wronskian yields I<sub>ν</sub>(x) directly. 1181 </p> 1182</div> 1183<table xmlns:rev="http://www.cs.rpi.edu/~gregod/boost/tools/doc/revision" width="100%"><tr> 1184<td align="left"></td> 1185<td align="right"><div class="copyright-footer">Copyright © 2006-2019 Nikhar 1186 Agrawal, Anton Bikineev, Paul A. Bristow, Marco Guazzone, Christopher Kormanyos, 1187 Hubert Holin, Bruno Lalande, John Maddock, Jeremy Murphy, Matthew Pulver, Johan 1188 Råde, Gautam Sewani, Benjamin Sobotta, Nicholas Thompson, Thijs van den Berg, 1189 Daryle Walker and Xiaogang Zhang<p> 1190 Distributed under the Boost Software License, Version 1.0. (See accompanying 1191 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>) 1192 </p> 1193</div></td> 1194</tr></table> 1195<hr> 1196<div class="spirit-nav"> 1197<a accesskey="p" href="bessel_root.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="sph_bessel.html"><img src="../../../../../../doc/src/images/next.png" alt="Next"></a> 1198</div> 1199</body> 1200</html> 1201
