1<html> 2<head> 3<meta http-equiv="Content-Type" content="text/html; charset=UTF-8"> 4<title>Finding Zeros of 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_first.html" title="Bessel Functions of the First and Second Kinds"> 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_first.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"> 26<div class="titlepage"><div><div><h3 class="title"> 27<a name="math_toolkit.bessel.bessel_root"></a><a class="link" href="bessel_root.html" title="Finding Zeros of Bessel Functions of the First and Second Kinds">Finding Zeros of Bessel 28 Functions of the First and Second Kinds</a> 29</h3></div></div></div> 30<h5> 31<a name="math_toolkit.bessel.bessel_root.h0"></a> 32 <span class="phrase"><a name="math_toolkit.bessel.bessel_root.synopsis"></a></span><a class="link" href="bessel_root.html#math_toolkit.bessel.bessel_root.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<p> 38 Functions for obtaining both a single zero or root of the Bessel function, 39 and placing multiple zeros into a container like <code class="computeroutput"><span class="identifier">std</span><span class="special">::</span><span class="identifier">vector</span></code> 40 by providing an output iterator. 41 </p> 42<p> 43 The signature of the single value functions are: 44 </p> 45<pre class="programlisting"><span class="keyword">template</span> <span class="special"><</span><span class="keyword">class</span> <span class="identifier">T</span><span class="special">></span> 46<span class="identifier">T</span> <span class="identifier">cyl_bessel_j_zero</span><span class="special">(</span> 47 <span class="identifier">T</span> <span class="identifier">v</span><span class="special">,</span> <span class="comment">// Floating-point value for Jv.</span> 48 <span class="keyword">int</span> <span class="identifier">m</span><span class="special">);</span> <span class="comment">// 1-based index of zero.</span> 49 50<span class="keyword">template</span> <span class="special"><</span><span class="keyword">class</span> <span class="identifier">T</span><span class="special">></span> 51<span class="identifier">T</span> <span class="identifier">cyl_neumann_zero</span><span class="special">(</span> 52 <span class="identifier">T</span> <span class="identifier">v</span><span class="special">,</span> <span class="comment">// Floating-point value for Jv.</span> 53 <span class="keyword">int</span> <span class="identifier">m</span><span class="special">);</span> <span class="comment">// 1-based index of zero.</span> 54</pre> 55<p> 56 and for multiple zeros: 57 </p> 58<pre class="programlisting"><span class="keyword">template</span> <span class="special"><</span><span class="keyword">class</span> <span class="identifier">T</span><span class="special">,</span> <span class="keyword">class</span> <span class="identifier">OutputIterator</span><span class="special">></span> 59<span class="identifier">OutputIterator</span> <span class="identifier">cyl_bessel_j_zero</span><span class="special">(</span> 60 <span class="identifier">T</span> <span class="identifier">v</span><span class="special">,</span> <span class="comment">// Floating-point value for Jv.</span> 61 <span class="keyword">int</span> <span class="identifier">start_index</span><span class="special">,</span> <span class="comment">// 1-based index of first zero.</span> 62 <span class="keyword">unsigned</span> <span class="identifier">number_of_zeros</span><span class="special">,</span> <span class="comment">// How many zeros to generate.</span> 63 <span class="identifier">OutputIterator</span> <span class="identifier">out_it</span><span class="special">);</span> <span class="comment">// Destination for zeros.</span> 64 65<span class="keyword">template</span> <span class="special"><</span><span class="keyword">class</span> <span class="identifier">T</span><span class="special">,</span> <span class="keyword">class</span> <span class="identifier">OutputIterator</span><span class="special">></span> 66<span class="identifier">OutputIterator</span> <span class="identifier">cyl_neumann_zero</span><span class="special">(</span> 67 <span class="identifier">T</span> <span class="identifier">v</span><span class="special">,</span> <span class="comment">// Floating-point value for Jv.</span> 68 <span class="keyword">int</span> <span class="identifier">start_index</span><span class="special">,</span> <span class="comment">// 1-based index of zero.</span> 69 <span class="keyword">unsigned</span> <span class="identifier">number_of_zeros</span><span class="special">,</span> <span class="comment">// How many zeros to generate</span> 70 <span class="identifier">OutputIterator</span> <span class="identifier">out_it</span><span class="special">);</span> <span class="comment">// Destination for zeros.</span> 71</pre> 72<p> 73 There are also versions which allow control of the <a class="link" href="../../policy.html" title="Chapter 21. Policies: Controlling Precision, Error Handling etc">Policies</a> 74 for error handling and precision. 75 </p> 76<pre class="programlisting"> <span class="keyword">template</span> <span class="special"><</span><span class="keyword">class</span> <span class="identifier">T</span><span class="special">></span> 77 <span class="identifier">T</span> <span class="identifier">cyl_bessel_j_zero</span><span class="special">(</span> 78 <span class="identifier">T</span> <span class="identifier">v</span><span class="special">,</span> <span class="comment">// Floating-point value for Jv.</span> 79 <span class="keyword">int</span> <span class="identifier">m</span><span class="special">,</span> <span class="comment">// 1-based index of zero.</span> 80 <span class="keyword">const</span> <span class="identifier">Policy</span><span class="special">&);</span> <span class="comment">// Policy to use.</span> 81 82 <span class="keyword">template</span> <span class="special"><</span><span class="keyword">class</span> <span class="identifier">T</span><span class="special">></span> 83 <span class="identifier">T</span> <span class="identifier">cyl_neumann_zero</span><span class="special">(</span> 84 <span class="identifier">T</span> <span class="identifier">v</span><span class="special">,</span> <span class="comment">// Floating-point value for Jv.</span> 85 <span class="keyword">int</span> <span class="identifier">m</span><span class="special">,</span> <span class="comment">// 1-based index of zero.</span> 86 <span class="keyword">const</span> <span class="identifier">Policy</span><span class="special">&);</span> <span class="comment">// Policy to use.</span> 87 88 89<span class="keyword">template</span> <span class="special"><</span><span class="keyword">class</span> <span class="identifier">T</span><span class="special">,</span> <span class="keyword">class</span> <span class="identifier">OutputIterator</span><span class="special">></span> 90<span class="identifier">OutputIterator</span> <span class="identifier">cyl_bessel_j_zero</span><span class="special">(</span> 91 <span class="identifier">T</span> <span class="identifier">v</span><span class="special">,</span> <span class="comment">// Floating-point value for Jv.</span> 92 <span class="keyword">int</span> <span class="identifier">start_index</span><span class="special">,</span> <span class="comment">// 1-based index of first zero.</span> 93 <span class="keyword">unsigned</span> <span class="identifier">number_of_zeros</span><span class="special">,</span> <span class="comment">// How many zeros to generate.</span> 94 <span class="identifier">OutputIterator</span> <span class="identifier">out_it</span><span class="special">,</span> <span class="comment">// Destination for zeros.</span> 95 <span class="keyword">const</span> <span class="identifier">Policy</span><span class="special">&</span> <span class="identifier">pol</span><span class="special">);</span> <span class="comment">// Policy to use.</span> 96 97<span class="keyword">template</span> <span class="special"><</span><span class="keyword">class</span> <span class="identifier">T</span><span class="special">,</span> <span class="keyword">class</span> <span class="identifier">OutputIterator</span><span class="special">></span> 98<span class="identifier">OutputIterator</span> <span class="identifier">cyl_neumann_zero</span><span class="special">(</span> 99 <span class="identifier">T</span> <span class="identifier">v</span><span class="special">,</span> <span class="comment">// Floating-point value for Jv.</span> 100 <span class="keyword">int</span> <span class="identifier">start_index</span><span class="special">,</span> <span class="comment">// 1-based index of zero.</span> 101 <span class="keyword">unsigned</span> <span class="identifier">number_of_zeros</span><span class="special">,</span> <span class="comment">// How many zeros to generate.</span> 102 <span class="identifier">OutputIterator</span> <span class="identifier">out_it</span><span class="special">,</span> <span class="comment">// Destination for zeros.</span> 103 <span class="keyword">const</span> <span class="identifier">Policy</span><span class="special">&</span> <span class="identifier">pol</span><span class="special">);</span> <span class="comment">// Policy to use.</span> 104</pre> 105<h5> 106<a name="math_toolkit.bessel.bessel_root.h1"></a> 107 <span class="phrase"><a name="math_toolkit.bessel.bessel_root.description"></a></span><a class="link" href="bessel_root.html#math_toolkit.bessel.bessel_root.description">Description</a> 108 </h5> 109<p> 110 Every real order ν cylindrical Bessel and Neumann functions have an infinite 111 number of zeros on the positive real axis. The real zeros on the positive 112 real axis can be found by solving for the roots of 113 </p> 114<div class="blockquote"><blockquote class="blockquote"><p> 115 <span class="emphasis"><em>J<sub>ν</sub>(j<sub>ν, m</sub>) = 0</em></span> 116 </p></blockquote></div> 117<div class="blockquote"><blockquote class="blockquote"><p> 118 <span class="emphasis"><em>Y<sub>ν</sub>(y<sub>ν, m</sub>) = 0</em></span> 119 </p></blockquote></div> 120<p> 121 Here, <span class="emphasis"><em>j<sub>ν, m</sub></em></span> represents the <span class="emphasis"><em>m<sup>th</sup></em></span> root 122 of the cylindrical Bessel function of order <span class="emphasis"><em>ν</em></span>, and <span class="emphasis"><em>y<sub>ν, 123 m</sub></em></span> represents the <span class="emphasis"><em>m<sup>th</sup></em></span> root of the cylindrical 124 Neumann function of order <span class="emphasis"><em>ν</em></span>. 125 </p> 126<p> 127 The zeros or roots (values of <code class="computeroutput"><span class="identifier">x</span></code> 128 where the function crosses the horizontal <code class="computeroutput"><span class="identifier">y</span> 129 <span class="special">=</span> <span class="number">0</span></code> 130 axis) of the Bessel and Neumann functions are computed by two functions, 131 <code class="computeroutput"><span class="identifier">cyl_bessel_j_zero</span></code> and <code class="computeroutput"><span class="identifier">cyl_neumann_zero</span></code>. 132 </p> 133<p> 134 In each case the index or rank of the zero returned is 1-based, which is 135 to say: 136 </p> 137<div class="blockquote"><blockquote class="blockquote"><p> 138 cyl_bessel_j_zero(v, 1); 139 </p></blockquote></div> 140<p> 141 returns the first zero of Bessel J. 142 </p> 143<p> 144 Passing an <code class="computeroutput"><span class="identifier">start_index</span> <span class="special"><=</span> 145 <span class="number">0</span></code> results in a <code class="computeroutput"><span class="identifier">std</span><span class="special">::</span><span class="identifier">domain_error</span></code> 146 being raised. 147 </p> 148<p> 149 For certain parameters, however, the zero'th root is defined and it has a 150 value of zero. For example, the zero'th root of <code class="computeroutput"><span class="identifier">J</span><span class="special">[</span><span class="identifier">sub</span> <span class="identifier">v</span><span class="special">](</span><span class="identifier">x</span><span class="special">)</span></code> 151 is defined and it has a value of zero for all values of <code class="computeroutput"><span class="identifier">v</span> 152 <span class="special">></span> <span class="number">0</span></code> 153 and for negative integer values of <code class="computeroutput"><span class="identifier">v</span> 154 <span class="special">=</span> <span class="special">-</span><span class="identifier">n</span></code>. Similar cases are described in the implementation 155 details below. 156 </p> 157<p> 158 The order <code class="computeroutput"><span class="identifier">v</span></code> of <code class="computeroutput"><span class="identifier">J</span></code> can be positive, negative and zero for 159 the <code class="computeroutput"><span class="identifier">cyl_bessel_j</span></code> and <code class="computeroutput"><span class="identifier">cyl_neumann</span></code> functions, but not infinite 160 nor NaN. 161 </p> 162<div class="blockquote"><blockquote class="blockquote"><p> 163 <span class="inlinemediaobject"><img src="../../../graphs/bessel_j_zeros.svg" align="middle"></span> 164 165 </p></blockquote></div> 166<div class="blockquote"><blockquote class="blockquote"><p> 167 <span class="inlinemediaobject"><img src="../../../graphs/neumann_y_zeros.svg" align="middle"></span> 168 169 </p></blockquote></div> 170<h5> 171<a name="math_toolkit.bessel.bessel_root.h2"></a> 172 <span class="phrase"><a name="math_toolkit.bessel.bessel_root.examples_of_finding_bessel_and_n"></a></span><a class="link" href="bessel_root.html#math_toolkit.bessel.bessel_root.examples_of_finding_bessel_and_n">Examples 173 of finding Bessel and Neumann zeros</a> 174 </h5> 175<p> 176 This example demonstrates calculating zeros of the Bessel and Neumann functions. 177 It also shows how Boost.Math and Boost.Multiprecision can be combined to 178 provide a many decimal digit precision. For 50 decimal digit precision we 179 need to include 180 </p> 181<pre class="programlisting"><span class="preprocessor">#include</span> <span class="special"><</span><span class="identifier">boost</span><span class="special">/</span><span class="identifier">multiprecision</span><span class="special">/</span><span class="identifier">cpp_dec_float</span><span class="special">.</span><span class="identifier">hpp</span><span class="special">></span> 182</pre> 183<p> 184 and a <code class="computeroutput"><span class="keyword">typedef</span></code> for <code class="computeroutput"><span class="identifier">float_type</span></code> may be convenient (allowing 185 a quick switch to re-compute at built-in <code class="computeroutput"><span class="keyword">double</span></code> 186 or other precision) 187 </p> 188<pre class="programlisting"><span class="keyword">typedef</span> <span class="identifier">boost</span><span class="special">::</span><span class="identifier">multiprecision</span><span class="special">::</span><span class="identifier">cpp_dec_float_50</span> <span class="identifier">float_type</span><span class="special">;</span> 189</pre> 190<p> 191 To use the functions for finding zeros of the functions we need 192 </p> 193<pre class="programlisting"><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> 194</pre> 195<p> 196 This file includes the forward declaration signatures for the zero-finding 197 functions: 198 </p> 199<pre class="programlisting"><span class="comment">// #include <boost/math/special_functions/math_fwd.hpp></span> 200</pre> 201<p> 202 but more details are in the full documentation, for example at <a href="http://www.boost.org/doc/libs/1_53_0/libs/math/doc/sf_and_dist/html/math_toolkit/special/bessel/bessel_over.html" target="_top">Boost.Math 203 Bessel functions</a>. 204 </p> 205<p> 206 This example shows obtaining both a single zero of the Bessel function, and 207 then placing multiple zeros into a container like <code class="computeroutput"><span class="identifier">std</span><span class="special">::</span><span class="identifier">vector</span></code> 208 by providing an iterator. 209 </p> 210<div class="tip"><table border="0" summary="Tip"> 211<tr> 212<td rowspan="2" align="center" valign="top" width="25"><img alt="[Tip]" src="../../../../../../doc/src/images/tip.png"></td> 213<th align="left">Tip</th> 214</tr> 215<tr><td align="left" valign="top"><p> 216 It is always wise to place code using Boost.Math inside try'n'catch blocks; 217 this will ensure that helpful error messages are shown when exceptional 218 conditions arise. 219 </p></td></tr> 220</table></div> 221<p> 222 First, evaluate a single Bessel zero. 223 </p> 224<p> 225 The precision is controlled by the float-point type of template parameter 226 <code class="computeroutput"><span class="identifier">T</span></code> of <code class="computeroutput"><span class="identifier">v</span></code> 227 so this example has <code class="computeroutput"><span class="keyword">double</span></code> precision, 228 at least 15 but up to 17 decimal digits (for the common 64-bit double). 229 </p> 230<pre class="programlisting"><span class="comment">// double root = boost::math::cyl_bessel_j_zero(0.0, 1);</span> 231<span class="comment">// // Displaying with default precision of 6 decimal digits:</span> 232<span class="comment">// std::cout << "boost::math::cyl_bessel_j_zero(0.0, 1) " << root << std::endl; // 2.40483</span> 233<span class="comment">// // And with all the guaranteed (15) digits:</span> 234<span class="comment">// std::cout.precision(std::numeric_limits<double>::digits10);</span> 235<span class="comment">// std::cout << "boost::math::cyl_bessel_j_zero(0.0, 1) " << root << std::endl; // 2.40482555769577</span> 236</pre> 237<p> 238 But note that because the parameter <code class="computeroutput"><span class="identifier">v</span></code> 239 controls the precision of the result, <code class="computeroutput"><span class="identifier">v</span></code> 240 <span class="bold"><strong>must be a floating-point type</strong></span>. So if you 241 provide an integer type, say 0, rather than 0.0, then it will fail to compile 242 thus: 243 </p> 244<pre class="programlisting"><span class="identifier">root</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">cyl_bessel_j_zero</span><span class="special">(</span><span class="number">0</span><span class="special">,</span> <span class="number">1</span><span class="special">);</span> 245</pre> 246<p> 247 with this error message 248 </p> 249<pre class="programlisting"><span class="identifier">error</span> <span class="identifier">C2338</span><span class="special">:</span> <span class="identifier">Order</span> <span class="identifier">must</span> <span class="identifier">be</span> <span class="identifier">a</span> <span class="identifier">floating</span><span class="special">-</span><span class="identifier">point</span> <span class="identifier">type</span><span class="special">.</span> 250</pre> 251<p> 252 Optionally, we can use a policy to ignore errors, C-style, returning some 253 value, perhaps infinity or NaN, or the best that can be done. (See <a class="link" href="../pol_tutorial/user_def_err_pol.html" title="Calling User Defined Error Handlers">user error handling</a>). 254 </p> 255<p> 256 To create a (possibly unwise!) policy <code class="computeroutput"><span class="identifier">ignore_all_policy</span></code> 257 that ignores all errors: 258 </p> 259<pre class="programlisting"><span class="keyword">typedef</span> <span class="identifier">boost</span><span class="special">::</span><span class="identifier">math</span><span class="special">::</span><span class="identifier">policies</span><span class="special">::</span><span class="identifier">policy</span><span class="special"><</span> 260 <span class="identifier">boost</span><span class="special">::</span><span class="identifier">math</span><span class="special">::</span><span class="identifier">policies</span><span class="special">::</span><span class="identifier">domain_error</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">policies</span><span class="special">::</span><span class="identifier">ignore_error</span><span class="special">>,</span> 261 <span class="identifier">boost</span><span class="special">::</span><span class="identifier">math</span><span class="special">::</span><span class="identifier">policies</span><span class="special">::</span><span class="identifier">overflow_error</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">policies</span><span class="special">::</span><span class="identifier">ignore_error</span><span class="special">>,</span> 262 <span class="identifier">boost</span><span class="special">::</span><span class="identifier">math</span><span class="special">::</span><span class="identifier">policies</span><span class="special">::</span><span class="identifier">underflow_error</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">policies</span><span class="special">::</span><span class="identifier">ignore_error</span><span class="special">>,</span> 263 <span class="identifier">boost</span><span class="special">::</span><span class="identifier">math</span><span class="special">::</span><span class="identifier">policies</span><span class="special">::</span><span class="identifier">denorm_error</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">policies</span><span class="special">::</span><span class="identifier">ignore_error</span><span class="special">>,</span> 264 <span class="identifier">boost</span><span class="special">::</span><span class="identifier">math</span><span class="special">::</span><span class="identifier">policies</span><span class="special">::</span><span class="identifier">pole_error</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">policies</span><span class="special">::</span><span class="identifier">ignore_error</span><span class="special">>,</span> 265 <span class="identifier">boost</span><span class="special">::</span><span class="identifier">math</span><span class="special">::</span><span class="identifier">policies</span><span class="special">::</span><span class="identifier">evaluation_error</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">policies</span><span class="special">::</span><span class="identifier">ignore_error</span><span class="special">></span> 266 <span class="special">></span> <span class="identifier">ignore_all_policy</span><span class="special">;</span> 267</pre> 268<p> 269 Examples of use of this <code class="computeroutput"><span class="identifier">ignore_all_policy</span></code> 270 are 271 </p> 272<pre class="programlisting"><span class="keyword">double</span> <span class="identifier">inf</span> <span class="special">=</span> <span class="identifier">std</span><span class="special">::</span><span class="identifier">numeric_limits</span><span class="special"><</span><span class="keyword">double</span><span class="special">>::</span><span class="identifier">infinity</span><span class="special">();</span> 273<span class="keyword">double</span> <span class="identifier">nan</span> <span class="special">=</span> <span class="identifier">std</span><span class="special">::</span><span class="identifier">numeric_limits</span><span class="special"><</span><span class="keyword">double</span><span class="special">>::</span><span class="identifier">quiet_NaN</span><span class="special">();</span> 274 275<span class="keyword">double</span> <span class="identifier">dodgy_root</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">cyl_bessel_j_zero</span><span class="special">(-</span><span class="number">1.0</span><span class="special">,</span> <span class="number">1</span><span class="special">,</span> <span class="identifier">ignore_all_policy</span><span class="special">());</span> 276<span class="identifier">std</span><span class="special">::</span><span class="identifier">cout</span> <span class="special"><<</span> <span class="string">"boost::math::cyl_bessel_j_zero(-1.0, 1) "</span> <span class="special"><<</span> <span class="identifier">dodgy_root</span> <span class="special"><<</span> <span class="identifier">std</span><span class="special">::</span><span class="identifier">endl</span><span class="special">;</span> <span class="comment">// 1.#QNAN</span> 277<span class="keyword">double</span> <span class="identifier">inf_root</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">cyl_bessel_j_zero</span><span class="special">(</span><span class="identifier">inf</span><span class="special">,</span> <span class="number">1</span><span class="special">,</span> <span class="identifier">ignore_all_policy</span><span class="special">());</span> 278<span class="identifier">std</span><span class="special">::</span><span class="identifier">cout</span> <span class="special"><<</span> <span class="string">"boost::math::cyl_bessel_j_zero(inf, 1) "</span> <span class="special"><<</span> <span class="identifier">inf_root</span> <span class="special"><<</span> <span class="identifier">std</span><span class="special">::</span><span class="identifier">endl</span><span class="special">;</span> <span class="comment">// 1.#QNAN</span> 279<span class="keyword">double</span> <span class="identifier">nan_root</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">cyl_bessel_j_zero</span><span class="special">(</span><span class="identifier">nan</span><span class="special">,</span> <span class="number">1</span><span class="special">,</span> <span class="identifier">ignore_all_policy</span><span class="special">());</span> 280<span class="identifier">std</span><span class="special">::</span><span class="identifier">cout</span> <span class="special"><<</span> <span class="string">"boost::math::cyl_bessel_j_zero(nan, 1) "</span> <span class="special"><<</span> <span class="identifier">nan_root</span> <span class="special"><<</span> <span class="identifier">std</span><span class="special">::</span><span class="identifier">endl</span><span class="special">;</span> <span class="comment">// 1.#QNAN</span> 281</pre> 282<p> 283 Another version of <code class="computeroutput"><span class="identifier">cyl_bessel_j_zero</span></code> 284 allows calculation of multiple zeros with one call, placing the results in 285 a container, often <code class="computeroutput"><span class="identifier">std</span><span class="special">::</span><span class="identifier">vector</span></code>. For example, generate and display 286 the first five <code class="computeroutput"><span class="keyword">double</span></code> roots 287 of J<sub>v</sub> for integral order 2, as column <span class="emphasis"><em>J<sub>2</sub>(x)</em></span> in table 288 1 of <a href="http://mathworld.wolfram.com/BesselFunctionZeros.html" target="_top">Wolfram 289 Bessel Function Zeros</a>. 290 </p> 291<pre class="programlisting"><span class="keyword">unsigned</span> <span class="keyword">int</span> <span class="identifier">n_roots</span> <span class="special">=</span> <span class="number">5U</span><span class="special">;</span> 292<span class="identifier">std</span><span class="special">::</span><span class="identifier">vector</span><span class="special"><</span><span class="keyword">double</span><span class="special">></span> <span class="identifier">roots</span><span class="special">;</span> 293<span class="identifier">boost</span><span class="special">::</span><span class="identifier">math</span><span class="special">::</span><span class="identifier">cyl_bessel_j_zero</span><span class="special">(</span><span class="number">2.0</span><span class="special">,</span> <span class="number">1</span><span class="special">,</span> <span class="identifier">n_roots</span><span class="special">,</span> <span class="identifier">std</span><span class="special">::</span><span class="identifier">back_inserter</span><span class="special">(</span><span class="identifier">roots</span><span class="special">));</span> 294<span class="identifier">std</span><span class="special">::</span><span class="identifier">copy</span><span class="special">(</span><span class="identifier">roots</span><span class="special">.</span><span class="identifier">begin</span><span class="special">(),</span> 295 <span class="identifier">roots</span><span class="special">.</span><span class="identifier">end</span><span class="special">(),</span> 296 <span class="identifier">std</span><span class="special">::</span><span class="identifier">ostream_iterator</span><span class="special"><</span><span class="keyword">double</span><span class="special">>(</span><span class="identifier">std</span><span class="special">::</span><span class="identifier">cout</span><span class="special">,</span> <span class="string">"\n"</span><span class="special">));</span> 297</pre> 298<p> 299 Or we can use Boost.Multiprecision to generate 50 decimal digit roots of 300 <span class="emphasis"><em>J<sub>v</sub></em></span> for non-integral order <code class="computeroutput"><span class="identifier">v</span><span class="special">=</span> <span class="number">71</span><span class="special">/</span><span class="number">19</span> <span class="special">==</span> <span class="number">3.736842</span></code>, 301 expressed as an exact-integer fraction to generate the most accurate value 302 possible for all floating-point types. 303 </p> 304<p> 305 We set the precision of the output stream, and show trailing zeros to display 306 a fixed 50 decimal digits. 307 </p> 308<pre class="programlisting"><span class="identifier">std</span><span class="special">::</span><span class="identifier">cout</span><span class="special">.</span><span class="identifier">precision</span><span class="special">(</span><span class="identifier">std</span><span class="special">::</span><span class="identifier">numeric_limits</span><span class="special"><</span><span class="identifier">float_type</span><span class="special">>::</span><span class="identifier">digits10</span><span class="special">);</span> <span class="comment">// 50 decimal digits.</span> 309<span class="identifier">std</span><span class="special">::</span><span class="identifier">cout</span> <span class="special"><<</span> <span class="identifier">std</span><span class="special">::</span><span class="identifier">showpoint</span> <span class="special"><<</span> <span class="identifier">std</span><span class="special">::</span><span class="identifier">endl</span><span class="special">;</span> <span class="comment">// Show trailing zeros.</span> 310 311<span class="identifier">float_type</span> <span class="identifier">x</span> <span class="special">=</span> <span class="identifier">float_type</span><span class="special">(</span><span class="number">71</span><span class="special">)</span> <span class="special">/</span> <span class="number">19</span><span class="special">;</span> 312<span class="identifier">float_type</span> <span class="identifier">r</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">cyl_bessel_j_zero</span><span class="special">(</span><span class="identifier">x</span><span class="special">,</span> <span class="number">1</span><span class="special">);</span> <span class="comment">// 1st root.</span> 313<span class="identifier">std</span><span class="special">::</span><span class="identifier">cout</span> <span class="special"><<</span> <span class="string">"x = "</span> <span class="special"><<</span> <span class="identifier">x</span> <span class="special"><<</span> <span class="string">", r = "</span> <span class="special"><<</span> <span class="identifier">r</span> <span class="special"><<</span> <span class="identifier">std</span><span class="special">::</span><span class="identifier">endl</span><span class="special">;</span> 314 315<span class="identifier">r</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">cyl_bessel_j_zero</span><span class="special">(</span><span class="identifier">x</span><span class="special">,</span> <span class="number">20U</span><span class="special">);</span> <span class="comment">// 20th root.</span> 316<span class="identifier">std</span><span class="special">::</span><span class="identifier">cout</span> <span class="special"><<</span> <span class="string">"x = "</span> <span class="special"><<</span> <span class="identifier">x</span> <span class="special"><<</span> <span class="string">", r = "</span> <span class="special"><<</span> <span class="identifier">r</span> <span class="special"><<</span> <span class="identifier">std</span><span class="special">::</span><span class="identifier">endl</span><span class="special">;</span> 317 318<span class="identifier">std</span><span class="special">::</span><span class="identifier">vector</span><span class="special"><</span><span class="identifier">float_type</span><span class="special">></span> <span class="identifier">zeros</span><span class="special">;</span> 319<span class="identifier">boost</span><span class="special">::</span><span class="identifier">math</span><span class="special">::</span><span class="identifier">cyl_bessel_j_zero</span><span class="special">(</span><span class="identifier">x</span><span class="special">,</span> <span class="number">1</span><span class="special">,</span> <span class="number">3</span><span class="special">,</span> <span class="identifier">std</span><span class="special">::</span><span class="identifier">back_inserter</span><span class="special">(</span><span class="identifier">zeros</span><span class="special">));</span> 320 321<span class="identifier">std</span><span class="special">::</span><span class="identifier">cout</span> <span class="special"><<</span> <span class="string">"cyl_bessel_j_zeros"</span> <span class="special"><<</span> <span class="identifier">std</span><span class="special">::</span><span class="identifier">endl</span><span class="special">;</span> 322<span class="comment">// Print the roots to the output stream.</span> 323<span class="identifier">std</span><span class="special">::</span><span class="identifier">copy</span><span class="special">(</span><span class="identifier">zeros</span><span class="special">.</span><span class="identifier">begin</span><span class="special">(),</span> <span class="identifier">zeros</span><span class="special">.</span><span class="identifier">end</span><span class="special">(),</span> 324 <span class="identifier">std</span><span class="special">::</span><span class="identifier">ostream_iterator</span><span class="special"><</span><span class="identifier">float_type</span><span class="special">>(</span><span class="identifier">std</span><span class="special">::</span><span class="identifier">cout</span><span class="special">,</span> <span class="string">"\n"</span><span class="special">));</span> 325</pre> 326<h6> 327<a name="math_toolkit.bessel.bessel_root.h3"></a> 328 <span class="phrase"><a name="math_toolkit.bessel.bessel_root.using_output_iterator_to_sum_zer"></a></span><a class="link" href="bessel_root.html#math_toolkit.bessel.bessel_root.using_output_iterator_to_sum_zer">Using 329 Output Iterator to sum zeros of Bessel Functions</a> 330 </h6> 331<p> 332 This example demonstrates summing zeros of the Bessel functions. To use the 333 functions for finding zeros of the functions we need 334 </p> 335<pre class="programlisting"><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> 336</pre> 337<p> 338 We use the <code class="computeroutput"><span class="identifier">cyl_bessel_j_zero</span></code> 339 output iterator parameter <code class="computeroutput"><span class="identifier">out_it</span></code> 340 to create a sum of <span class="emphasis"><em>1/zeros<sup>2</sup></em></span> by defining a custom output 341 iterator: 342 </p> 343<pre class="programlisting"><span class="keyword">template</span> <span class="special"><</span><span class="keyword">class</span> <span class="identifier">T</span><span class="special">></span> 344<span class="keyword">struct</span> <span class="identifier">output_summation_iterator</span> 345<span class="special">{</span> 346 <span class="identifier">output_summation_iterator</span><span class="special">(</span><span class="identifier">T</span><span class="special">*</span> <span class="identifier">p</span><span class="special">)</span> <span class="special">:</span> <span class="identifier">p_sum</span><span class="special">(</span><span class="identifier">p</span><span class="special">)</span> 347 <span class="special">{}</span> 348 <span class="identifier">output_summation_iterator</span><span class="special">&</span> <span class="keyword">operator</span><span class="special">*()</span> 349 <span class="special">{</span> <span class="keyword">return</span> <span class="special">*</span><span class="keyword">this</span><span class="special">;</span> <span class="special">}</span> 350 <span class="identifier">output_summation_iterator</span><span class="special">&</span> <span class="keyword">operator</span><span class="special">++()</span> 351 <span class="special">{</span> <span class="keyword">return</span> <span class="special">*</span><span class="keyword">this</span><span class="special">;</span> <span class="special">}</span> 352 <span class="identifier">output_summation_iterator</span><span class="special">&</span> <span class="keyword">operator</span><span class="special">++(</span><span class="keyword">int</span><span class="special">)</span> 353 <span class="special">{</span> <span class="keyword">return</span> <span class="special">*</span><span class="keyword">this</span><span class="special">;</span> <span class="special">}</span> 354 <span class="identifier">output_summation_iterator</span><span class="special">&</span> <span class="keyword">operator</span> <span class="special">=</span> <span class="special">(</span><span class="identifier">T</span> <span class="keyword">const</span><span class="special">&</span> <span class="identifier">val</span><span class="special">)</span> 355 <span class="special">{</span> 356 <span class="special">*</span><span class="identifier">p_sum</span> <span class="special">+=</span> <span class="number">1.</span><span class="special">/</span> <span class="special">(</span><span class="identifier">val</span> <span class="special">*</span> <span class="identifier">val</span><span class="special">);</span> <span class="comment">// Summing 1/zero^2.</span> 357 <span class="keyword">return</span> <span class="special">*</span><span class="keyword">this</span><span class="special">;</span> 358 <span class="special">}</span> 359<span class="keyword">private</span><span class="special">:</span> 360 <span class="identifier">T</span><span class="special">*</span> <span class="identifier">p_sum</span><span class="special">;</span> 361<span class="special">};</span> 362</pre> 363<p> 364 The sum is calculated for many values, converging on the analytical exact 365 value of <code class="computeroutput"><span class="number">1</span><span class="special">/</span><span class="number">8</span></code>. 366 </p> 367<pre class="programlisting"><span class="keyword">using</span> <span class="identifier">boost</span><span class="special">::</span><span class="identifier">math</span><span class="special">::</span><span class="identifier">cyl_bessel_j_zero</span><span class="special">;</span> 368<span class="keyword">double</span> <span class="identifier">nu</span> <span class="special">=</span> <span class="number">1.</span><span class="special">;</span> 369<span class="keyword">double</span> <span class="identifier">sum</span> <span class="special">=</span> <span class="number">0</span><span class="special">;</span> 370<span class="identifier">output_summation_iterator</span><span class="special"><</span><span class="keyword">double</span><span class="special">></span> <span class="identifier">it</span><span class="special">(&</span><span class="identifier">sum</span><span class="special">);</span> <span class="comment">// sum of 1/zeros^2</span> 371<span class="identifier">cyl_bessel_j_zero</span><span class="special">(</span><span class="identifier">nu</span><span class="special">,</span> <span class="number">1</span><span class="special">,</span> <span class="number">10000</span><span class="special">,</span> <span class="identifier">it</span><span class="special">);</span> 372 373<span class="keyword">double</span> <span class="identifier">s</span> <span class="special">=</span> <span class="number">1</span><span class="special">/(</span><span class="number">4</span> <span class="special">*</span> <span class="special">(</span><span class="identifier">nu</span> <span class="special">+</span> <span class="number">1</span><span class="special">));</span> <span class="comment">// 0.125 = 1/8 is exact analytical solution.</span> 374<span class="identifier">std</span><span class="special">::</span><span class="identifier">cout</span> <span class="special"><<</span> <span class="identifier">std</span><span class="special">::</span><span class="identifier">setprecision</span><span class="special">(</span><span class="number">6</span><span class="special">)</span> <span class="special"><<</span> <span class="string">"nu = "</span> <span class="special"><<</span> <span class="identifier">nu</span> <span class="special"><<</span> <span class="string">", sum = "</span> <span class="special"><<</span> <span class="identifier">sum</span> 375 <span class="special"><<</span> <span class="string">", exact = "</span> <span class="special"><<</span> <span class="identifier">s</span> <span class="special"><<</span> <span class="identifier">std</span><span class="special">::</span><span class="identifier">endl</span><span class="special">;</span> 376<span class="comment">// nu = 1.00000, sum = 0.124990, exact = 0.125000</span> 377</pre> 378<h6> 379<a name="math_toolkit.bessel.bessel_root.h4"></a> 380 <span class="phrase"><a name="math_toolkit.bessel.bessel_root.calculating_zeros_of_the_neumann"></a></span><a class="link" href="bessel_root.html#math_toolkit.bessel.bessel_root.calculating_zeros_of_the_neumann">Calculating 381 zeros of the Neumann function.</a> 382 </h6> 383<p> 384 This example also shows how Boost.Math and Boost.Multiprecision can be combined 385 to provide a many decimal digit precision. For 50 decimal digit precision 386 we need to include 387 </p> 388<pre class="programlisting"><span class="preprocessor">#include</span> <span class="special"><</span><span class="identifier">boost</span><span class="special">/</span><span class="identifier">multiprecision</span><span class="special">/</span><span class="identifier">cpp_dec_float</span><span class="special">.</span><span class="identifier">hpp</span><span class="special">></span> 389</pre> 390<p> 391 and a <code class="computeroutput"><span class="keyword">typedef</span></code> for <code class="computeroutput"><span class="identifier">float_type</span></code> may be convenient (allowing 392 a quick switch to re-compute at built-in <code class="computeroutput"><span class="keyword">double</span></code> 393 or other precision) 394 </p> 395<pre class="programlisting"><span class="keyword">typedef</span> <span class="identifier">boost</span><span class="special">::</span><span class="identifier">multiprecision</span><span class="special">::</span><span class="identifier">cpp_dec_float_50</span> <span class="identifier">float_type</span><span class="special">;</span> 396</pre> 397<p> 398 To use the functions for finding zeros of the <code class="computeroutput"><span class="identifier">cyl_neumann</span></code> 399 function we need: 400 </p> 401<pre class="programlisting"><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> 402</pre> 403<p> 404 The Neumann (Bessel Y) function zeros are evaluated very similarly: 405 </p> 406<pre class="programlisting"><span class="keyword">using</span> <span class="identifier">boost</span><span class="special">::</span><span class="identifier">math</span><span class="special">::</span><span class="identifier">cyl_neumann_zero</span><span class="special">;</span> 407<span class="keyword">double</span> <span class="identifier">zn</span> <span class="special">=</span> <span class="identifier">cyl_neumann_zero</span><span class="special">(</span><span class="number">2.</span><span class="special">,</span> <span class="number">1</span><span class="special">);</span> 408<span class="identifier">std</span><span class="special">::</span><span class="identifier">cout</span> <span class="special"><<</span> <span class="string">"cyl_neumann_zero(2., 1) = "</span> <span class="special"><<</span> <span class="identifier">zn</span> <span class="special"><<</span> <span class="identifier">std</span><span class="special">::</span><span class="identifier">endl</span><span class="special">;</span> 409 410<span class="identifier">std</span><span class="special">::</span><span class="identifier">vector</span><span class="special"><</span><span class="keyword">float</span><span class="special">></span> <span class="identifier">nzeros</span><span class="special">(</span><span class="number">3</span><span class="special">);</span> <span class="comment">// Space for 3 zeros.</span> 411<span class="identifier">cyl_neumann_zero</span><span class="special"><</span><span class="keyword">float</span><span class="special">>(</span><span class="number">2.F</span><span class="special">,</span> <span class="number">1</span><span class="special">,</span> <span class="identifier">nzeros</span><span class="special">.</span><span class="identifier">size</span><span class="special">(),</span> <span class="identifier">nzeros</span><span class="special">.</span><span class="identifier">begin</span><span class="special">());</span> 412 413<span class="identifier">std</span><span class="special">::</span><span class="identifier">cout</span> <span class="special"><<</span> <span class="string">"cyl_neumann_zero<float>(2.F, 1, "</span><span class="special">;</span> 414<span class="comment">// Print the zeros to the output stream.</span> 415<span class="identifier">std</span><span class="special">::</span><span class="identifier">copy</span><span class="special">(</span><span class="identifier">nzeros</span><span class="special">.</span><span class="identifier">begin</span><span class="special">(),</span> <span class="identifier">nzeros</span><span class="special">.</span><span class="identifier">end</span><span class="special">(),</span> 416 <span class="identifier">std</span><span class="special">::</span><span class="identifier">ostream_iterator</span><span class="special"><</span><span class="keyword">float</span><span class="special">>(</span><span class="identifier">std</span><span class="special">::</span><span class="identifier">cout</span><span class="special">,</span> <span class="string">", "</span><span class="special">));</span> 417 418<span class="identifier">std</span><span class="special">::</span><span class="identifier">cout</span> <span class="special"><<</span> <span class="string">"\n"</span><span class="string">"cyl_neumann_zero(static_cast<float_type>(220)/100, 1) = "</span> 419 <span class="special"><<</span> <span class="identifier">cyl_neumann_zero</span><span class="special">(</span><span class="keyword">static_cast</span><span class="special"><</span><span class="identifier">float_type</span><span class="special">>(</span><span class="number">220</span><span class="special">)/</span><span class="number">100</span><span class="special">,</span> <span class="number">1</span><span class="special">)</span> <span class="special"><<</span> <span class="identifier">std</span><span class="special">::</span><span class="identifier">endl</span><span class="special">;</span> 420<span class="comment">// 3.6154383428745996706772556069431792744372398748422</span> 421</pre> 422<h6> 423<a name="math_toolkit.bessel.bessel_root.h5"></a> 424 <span class="phrase"><a name="math_toolkit.bessel.bessel_root.error_messages_from_bad_input"></a></span><a class="link" href="bessel_root.html#math_toolkit.bessel.bessel_root.error_messages_from_bad_input">Error 425 messages from 'bad' input</a> 426 </h6> 427<p> 428 Another example demonstrates calculating zeros of the Bessel functions showing 429 the error messages from 'bad' input is handled by throwing exceptions. 430 </p> 431<p> 432 To use the functions for finding zeros of the functions we need: 433 </p> 434<pre class="programlisting"><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> 435<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">airy</span><span class="special">.</span><span class="identifier">hpp</span><span class="special">></span> 436</pre> 437<div class="tip"><table border="0" summary="Tip"> 438<tr> 439<td rowspan="2" align="center" valign="top" width="25"><img alt="[Tip]" src="../../../../../../doc/src/images/tip.png"></td> 440<th align="left">Tip</th> 441</tr> 442<tr><td align="left" valign="top"><p> 443 It is always wise to place all code using Boost.Math inside try'n'catch 444 blocks; this will ensure that helpful error messages can be shown when 445 exceptional conditions arise. 446 </p></td></tr> 447</table></div> 448<p> 449 Examples below show messages from several 'bad' arguments that throw a <code class="computeroutput"><span class="identifier">domain_error</span></code> exception. 450 </p> 451<pre class="programlisting"><span class="keyword">try</span> 452<span class="special">{</span> <span class="comment">// Try a zero order v.</span> 453 <span class="keyword">float</span> <span class="identifier">dodgy_root</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">cyl_bessel_j_zero</span><span class="special">(</span><span class="number">0.F</span><span class="special">,</span> <span class="number">0</span><span class="special">);</span> 454 <span class="identifier">std</span><span class="special">::</span><span class="identifier">cout</span> <span class="special"><<</span> <span class="string">"boost::math::cyl_bessel_j_zero(0.F, 0) "</span> <span class="special"><<</span> <span class="identifier">dodgy_root</span> <span class="special"><<</span> <span class="identifier">std</span><span class="special">::</span><span class="identifier">endl</span><span class="special">;</span> 455 <span class="comment">// Thrown exception Error in function boost::math::cyl_bessel_j_zero<double>(double, int):</span> 456 <span class="comment">// Requested the 0'th zero of J0, but the rank must be > 0 !</span> 457<span class="special">}</span> 458<span class="keyword">catch</span> <span class="special">(</span><span class="identifier">std</span><span class="special">::</span><span class="identifier">exception</span><span class="special">&</span> <span class="identifier">ex</span><span class="special">)</span> 459<span class="special">{</span> 460 <span class="identifier">std</span><span class="special">::</span><span class="identifier">cout</span> <span class="special"><<</span> <span class="string">"Thrown exception "</span> <span class="special"><<</span> <span class="identifier">ex</span><span class="special">.</span><span class="identifier">what</span><span class="special">()</span> <span class="special"><<</span> <span class="identifier">std</span><span class="special">::</span><span class="identifier">endl</span><span class="special">;</span> 461<span class="special">}</span> 462</pre> 463<div class="note"><table border="0" summary="Note"> 464<tr> 465<td rowspan="2" align="center" valign="top" width="25"><img alt="[Note]" src="../../../../../../doc/src/images/note.png"></td> 466<th align="left">Note</th> 467</tr> 468<tr><td align="left" valign="top"><p> 469 The type shown in the error message is the type <span class="bold"><strong>after 470 promotion</strong></span>, using <a class="link" href="../pol_ref/precision_pol.html" title="Precision Policies">precision 471 policy</a> and <a class="link" href="../pol_ref/internal_promotion.html" title="Internal Floating-point Promotion Policies">internal 472 promotion policy</a>, from <code class="computeroutput"><span class="keyword">float</span></code> 473 to <code class="computeroutput"><span class="keyword">double</span></code> in this case. 474 </p></td></tr> 475</table></div> 476<p> 477 In this example the promotion goes: 478 </p> 479<div class="orderedlist"><ol class="orderedlist" type="1"> 480<li class="listitem"> 481 Arguments are <code class="computeroutput"><span class="keyword">float</span></code> and 482 <code class="computeroutput"><span class="keyword">int</span></code>. 483 </li> 484<li class="listitem"> 485 Treat <code class="computeroutput"><span class="keyword">int</span></code> "as if" 486 it were a <code class="computeroutput"><span class="keyword">double</span></code>, so arguments 487 are <code class="computeroutput"><span class="keyword">float</span></code> and <code class="computeroutput"><span class="keyword">double</span></code>. 488 </li> 489<li class="listitem"> 490 Common type is <code class="computeroutput"><span class="keyword">double</span></code> - 491 so that's the precision we want (and the type that will be returned). 492 </li> 493<li class="listitem"> 494 Evaluate internally as <code class="computeroutput"><span class="keyword">double</span></code> 495 for full <code class="computeroutput"><span class="keyword">float</span></code> precision. 496 </li> 497</ol></div> 498<p> 499 See full code for other examples that promote from <code class="computeroutput"><span class="keyword">double</span></code> 500 to <code class="computeroutput"><span class="keyword">long</span> <span class="keyword">double</span></code>. 501 </p> 502<p> 503 Other examples of 'bad' inputs like infinity and NaN are below. Some compiler 504 warnings indicate that 'bad' values are detected at compile time. 505 </p> 506<pre class="programlisting"><span class="keyword">try</span> 507<span class="special">{</span> <span class="comment">// order v = inf</span> 508 <span class="identifier">std</span><span class="special">::</span><span class="identifier">cout</span> <span class="special"><<</span> <span class="string">"boost::math::cyl_bessel_j_zero(inf, 1) "</span> <span class="special"><<</span> <span class="identifier">std</span><span class="special">::</span><span class="identifier">endl</span><span class="special">;</span> 509 <span class="keyword">double</span> <span class="identifier">inf</span> <span class="special">=</span> <span class="identifier">std</span><span class="special">::</span><span class="identifier">numeric_limits</span><span class="special"><</span><span class="keyword">double</span><span class="special">>::</span><span class="identifier">infinity</span><span class="special">();</span> 510 <span class="keyword">double</span> <span class="identifier">inf_root</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">cyl_bessel_j_zero</span><span class="special">(</span><span class="identifier">inf</span><span class="special">,</span> <span class="number">1</span><span class="special">);</span> 511 <span class="identifier">std</span><span class="special">::</span><span class="identifier">cout</span> <span class="special"><<</span> <span class="string">"boost::math::cyl_bessel_j_zero(inf, 1) "</span> <span class="special"><<</span> <span class="identifier">inf_root</span> <span class="special"><<</span> <span class="identifier">std</span><span class="special">::</span><span class="identifier">endl</span><span class="special">;</span> 512 <span class="comment">// Throw exception Error in function boost::math::cyl_bessel_j_zero<long double>(long double, unsigned):</span> 513 <span class="comment">// Order argument is 1.#INF, but must be finite >= 0 !</span> 514<span class="special">}</span> 515<span class="keyword">catch</span> <span class="special">(</span><span class="identifier">std</span><span class="special">::</span><span class="identifier">exception</span><span class="special">&</span> <span class="identifier">ex</span><span class="special">)</span> 516<span class="special">{</span> 517 <span class="identifier">std</span><span class="special">::</span><span class="identifier">cout</span> <span class="special"><<</span> <span class="string">"Thrown exception "</span> <span class="special"><<</span> <span class="identifier">ex</span><span class="special">.</span><span class="identifier">what</span><span class="special">()</span> <span class="special"><<</span> <span class="identifier">std</span><span class="special">::</span><span class="identifier">endl</span><span class="special">;</span> 518<span class="special">}</span> 519 520<span class="keyword">try</span> 521<span class="special">{</span> <span class="comment">// order v = NaN, rank m = 1</span> 522 <span class="identifier">std</span><span class="special">::</span><span class="identifier">cout</span> <span class="special"><<</span> <span class="string">"boost::math::cyl_bessel_j_zero(nan, 1) "</span> <span class="special"><<</span> <span class="identifier">std</span><span class="special">::</span><span class="identifier">endl</span><span class="special">;</span> 523 <span class="keyword">double</span> <span class="identifier">nan</span> <span class="special">=</span> <span class="identifier">std</span><span class="special">::</span><span class="identifier">numeric_limits</span><span class="special"><</span><span class="keyword">double</span><span class="special">>::</span><span class="identifier">quiet_NaN</span><span class="special">();</span> 524 <span class="keyword">double</span> <span class="identifier">nan_root</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">cyl_bessel_j_zero</span><span class="special">(</span><span class="identifier">nan</span><span class="special">,</span> <span class="number">1</span><span class="special">);</span> 525 <span class="identifier">std</span><span class="special">::</span><span class="identifier">cout</span> <span class="special"><<</span> <span class="string">"boost::math::cyl_bessel_j_zero(nan, 1) "</span> <span class="special"><<</span> <span class="identifier">nan_root</span> <span class="special"><<</span> <span class="identifier">std</span><span class="special">::</span><span class="identifier">endl</span><span class="special">;</span> 526 <span class="comment">// Throw exception Error in function boost::math::cyl_bessel_j_zero<long double>(long double, unsigned):</span> 527 <span class="comment">// Order argument is 1.#QNAN, but must be finite >= 0 !</span> 528<span class="special">}</span> 529<span class="keyword">catch</span> <span class="special">(</span><span class="identifier">std</span><span class="special">::</span><span class="identifier">exception</span><span class="special">&</span> <span class="identifier">ex</span><span class="special">)</span> 530<span class="special">{</span> 531 <span class="identifier">std</span><span class="special">::</span><span class="identifier">cout</span> <span class="special"><<</span> <span class="string">"Thrown exception "</span> <span class="special"><<</span> <span class="identifier">ex</span><span class="special">.</span><span class="identifier">what</span><span class="special">()</span> <span class="special"><<</span> <span class="identifier">std</span><span class="special">::</span><span class="identifier">endl</span><span class="special">;</span> 532<span class="special">}</span> 533</pre> 534<p> 535 The output from other examples are shown appended to the full code listing. 536 </p> 537<p> 538 The full code (and output) for these examples is at <a href="../../../../example/bessel_zeros_example_1.cpp" target="_top">Bessel 539 zeros</a>, <a href="../../../../example/bessel_zeros_interator_example.cpp" target="_top">Bessel 540 zeros iterator</a>, <a href="../../../../example/neumann_zeros_example_1.cpp" target="_top">Neumann 541 zeros</a>, <a href="../../../../example/bessel_errors_example.cpp" target="_top">Bessel 542 error messages</a>. 543 </p> 544<h4> 545<a name="math_toolkit.bessel.bessel_root.h6"></a> 546 <span class="phrase"><a name="math_toolkit.bessel.bessel_root.implementation"></a></span><a class="link" href="bessel_root.html#math_toolkit.bessel.bessel_root.implementation">Implementation</a> 547 </h4> 548<p> 549 Various methods are used to compute initial estimates for <span class="emphasis"><em>j<sub>ν, m</sub></em></span> 550 and <span class="emphasis"><em>y<sub>ν, m</sub></em></span> ; these are described in detail below. 551 </p> 552<p> 553 After finding the initial estimate of a given root, its precision is subsequently 554 refined to the desired level using Newton-Raphson iteration from Boost.Math's 555 <a class="link" href="../roots_deriv.html" title="Root Finding With Derivatives: Newton-Raphson, Halley & Schröder">root-finding with derivatives</a> 556 utilities combined with the functions <a class="link" href="bessel_first.html" title="Bessel Functions of the First and Second Kinds">cyl_bessel_j</a> 557 and <a class="link" href="bessel_first.html" title="Bessel Functions of the First and Second Kinds">cyl_neumann</a>. 558 </p> 559<p> 560 Newton iteration requires both <span class="emphasis"><em>J<sub>ν</sub>(x)</em></span> or <span class="emphasis"><em>Y<sub>ν</sub>(x)</em></span> 561 as well as its derivative. The derivatives of <span class="emphasis"><em>J<sub>ν</sub>(x)</em></span> and 562 <span class="emphasis"><em>Y<sub>ν</sub>(x)</em></span> with respect to <span class="emphasis"><em>x</em></span> are given 563 by M. Abramowitz and I. A. Stegun, Handbook of Mathematical Functions, NBS 564 (1964). In particular, 565 </p> 566<div class="blockquote"><blockquote class="blockquote"><p> 567 <span class="serif_italic">d/<sub>dx</sub> <span class="emphasis"><em>J<sub>ν</sub>(x)</em></span> = <span class="emphasis"><em>J<sub>ν-1</sub>(x)</em></span> 568 - ν J<sub>ν</sub>(x) / x</span> 569 </p></blockquote></div> 570<div class="blockquote"><blockquote class="blockquote"><p> 571 <span class="serif_italic">d/<sub>dx</sub> <span class="emphasis"><em>Y<sub>ν</sub>(x)</em></span> = <span class="emphasis"><em>Y<sub>ν-1</sub>(x)</em></span> 572 - ν Y<sub>ν</sub>(x) / x</span> 573 </p></blockquote></div> 574<p> 575 Enumeration of the rank of a root (in other words the index of a root) begins 576 with one and counts up, in other words <span class="emphasis"><em>m,=1,2,3,…</em></span> The 577 value of the first root is always greater than zero. 578 </p> 579<p> 580 For certain special parameters, cylindrical Bessel functions and cylindrical 581 Neumann functions have a root at the origin. For example, <span class="emphasis"><em>J<sub>ν</sub>(x)</em></span> 582 has a root at the origin for every positive order <span class="emphasis"><em>ν > 0</em></span>, 583 and for every negative integer order <span class="emphasis"><em>ν = -n</em></span> with <span class="emphasis"><em>n 584 ∈ ℕ <sup>+</sup></em></span> and <span class="emphasis"><em>n ≠ 0</em></span>. 585 </p> 586<p> 587 In addition, <span class="emphasis"><em>Y<sub>ν</sub>(x)</em></span> has a root at the origin for every 588 negative half-integer order <span class="emphasis"><em>ν = -n/2</em></span>, with <span class="emphasis"><em>n 589 ∈ ℕ <sup>+</sup></em></span> and and <span class="emphasis"><em>n ≠ 0</em></span>. 590 </p> 591<p> 592 For these special parameter values, the origin with a value of <span class="emphasis"><em>x 593 = 0</em></span> is provided as the <span class="emphasis"><em>0<sup>th</sup></em></span> root generated 594 by <code class="computeroutput"><span class="identifier">cyl_bessel_j_zero</span><span class="special">()</span></code> 595 and <code class="computeroutput"><span class="identifier">cyl_neumann_zero</span><span class="special">()</span></code>. 596 </p> 597<p> 598 When calculating initial estimates for the roots of Bessel functions, a distinction 599 is made between positive order and negative order, and different methods 600 are used for these. In addition, different algorithms are used for the first 601 root <span class="emphasis"><em>m = 1</em></span> and for subsequent roots with higher rank 602 <span class="emphasis"><em>m ≥ 2</em></span>. Furthermore, estimates of the roots for Bessel 603 functions with order above and below a cutoff at <span class="emphasis"><em>ν = 2.2</em></span> 604 are calculated with different methods. 605 </p> 606<p> 607 Calculations of the estimates of <span class="emphasis"><em>j<sub>ν,1</sub></em></span> and <span class="emphasis"><em>y<sub>ν,1</sub></em></span> 608 with <span class="emphasis"><em>0 ≤ ν < 2.2</em></span> use empirically tabulated values. The 609 coefficients for these have been generated by a computer algebra system. 610 </p> 611<p> 612 Calculations of the estimates of <span class="emphasis"><em>j<sub>ν,1</sub></em></span> and <span class="emphasis"><em>y<sub>ν,1</sub></em></span> 613 with <span class="emphasis"><em>ν≥ 2.2</em></span> use Eqs.9.5.14 and 9.5.15 in M. Abramowitz 614 and I. A. Stegun, Handbook of Mathematical Functions, NBS (1964). 615 </p> 616<p> 617 In particular, 618 </p> 619<div class="blockquote"><blockquote class="blockquote"><p> 620 <span class="serif_italic">j<sub>ν,1</sub> ≅ ν + 1.85575 ν<sup>⅓</sup> + 1.033150 ν<sup>-⅓</sup> - 0.00397 ν<sup>-1</sup> - 0.0908 621 ν<sup>-5/3</sup> + 0.043 ν<sup>-7/3</sup> + …</span> 622 </p></blockquote></div> 623<p> 624 and 625 </p> 626<div class="blockquote"><blockquote class="blockquote"><p> 627 <span class="serif_italic">y<sub>ν,1</sub> ≅ ν + 0.93157 ν<sup>⅓</sup> + 0.26035 ν<sup>-⅓</sup> + 0.01198 ν<sup>-1</sup> - 0.0060 628 ν<sup>-5/3</sup> - 0.001 ν<sup>-7/3</sup> + …</span> 629 </p></blockquote></div> 630<p> 631 Calculations of the estimates of <span class="emphasis"><em>j<sub>ν, m</sub></em></span> and <span class="emphasis"><em>y<sub>ν, 632 m</sub></em></span> with rank <span class="emphasis"><em>m > 2</em></span> and <span class="emphasis"><em>0 ≤ ν < 633 2.2</em></span> use McMahon's approximation, as described in M. Abramowitz 634 and I. A. Stegan, Section 9.5 and 9.5.12. In particular, 635 </p> 636<div class="blockquote"><blockquote class="blockquote"><p> 637 <span class="emphasis"><em>j<sub>ν,m</sub>, y<sub>ν,m</sub> ≅</em></span> 638 </p></blockquote></div> 639<div class="blockquote"><blockquote class="blockquote"><div class="blockquote"><blockquote class="blockquote"><p> 640 β - (μ-1) / 8β 641 </p></blockquote></div></blockquote></div> 642<div class="blockquote"><blockquote class="blockquote"><div class="blockquote"><blockquote class="blockquote"><p> 643 <span class="emphasis"><em>- 4(μ-1)(7μ - 31) / 3(8β)<sup>3</sup></em></span> 644 </p></blockquote></div></blockquote></div> 645<div class="blockquote"><blockquote class="blockquote"><div class="blockquote"><blockquote class="blockquote"><p> 646 <span class="emphasis"><em>-32(μ-1)(83μ² - 982μ + 3779) / 15(8β)<sup>5</sup></em></span> 647 </p></blockquote></div></blockquote></div> 648<div class="blockquote"><blockquote class="blockquote"><div class="blockquote"><blockquote class="blockquote"><p> 649 <span class="emphasis"><em>-64(μ-1)(6949μ<sup>3</sup> - 153855μ² + 1585743μ- 6277237) / 105(8a)<sup>7</sup></em></span> 650 </p></blockquote></div></blockquote></div> 651<div class="blockquote"><blockquote class="blockquote"><div class="blockquote"><blockquote class="blockquote"><p> 652 <span class="emphasis"><em>- …</em></span> (5) 653 </p></blockquote></div></blockquote></div> 654<p> 655 where <span class="emphasis"><em>μ = 4ν<sup>2</sup></em></span> and <span class="emphasis"><em>β = (m + ½ν - ¼)π</em></span> for 656 <span class="emphasis"><em>j<sub>ν,m</sub></em></span> and <span class="emphasis"><em>β = (m + ½ν -¾)π for <span class="emphasis"><em>y<sub>ν,m</sub></em></span></em></span>. 657 </p> 658<p> 659 Calculations of the estimates of <span class="emphasis"><em>j<sub>ν, m</sub></em></span> and <span class="emphasis"><em>y<sub>ν, 660 m</sub></em></span> with <span class="emphasis"><em>ν ≥ 2.2</em></span> use one term in the asymptotic 661 expansion given in Eq.9.5.22 and top line of Eq.9.5.26 combined with Eq. 662 9.3.39, all in M. Abramowitz and I. A. Stegun, Handbook of Mathematical Functions, 663 NBS (1964) explicit and easy-to-understand treatment for asymptotic expansion 664 of zeros. The latter two equations are expressed for argument <span class="emphasis"><em>(x)</em></span> 665 greater than one. (Olver also gives the series form of the equations in 666 <a href="http://dlmf.nist.gov/10.21#vi" target="_top">§10.21(vi) McMahon's Asymptotic 667 Expansions for Large Zeros</a> - using slightly different variable names). 668 </p> 669<p> 670 In summary, 671 </p> 672<div class="blockquote"><blockquote class="blockquote"><p> 673 <span class="serif_italic">j<sub>ν, m</sub> ∼ νx(-ζ) + f<sub>1</sub>(-ζ/ν)</span> 674 </p></blockquote></div> 675<p> 676 where <span class="emphasis"><em>-ζ = ν<sup>-2/3</sup>a<sub>m</sub></em></span> and <span class="emphasis"><em>a<sub>m</sub></em></span> is the absolute 677 value of the <span class="emphasis"><em>m<sup>th</sup></em></span> root of <span class="emphasis"><em>Ai(x)</em></span> 678 on the negative real axis. 679 </p> 680<p> 681 Here <span class="emphasis"><em>x = x(-ζ)</em></span> is the inverse of the function 682 </p> 683<div class="blockquote"><blockquote class="blockquote"><p> 684 <span class="serif_italic">⅔(-ζ)<sup>3/2</sup> = √(x² - 1) - cos⁻¹(1/x)</span> 685 </p></blockquote></div> 686<p> 687 (7) 688 </p> 689<p> 690 Furthermore, 691 </p> 692<div class="blockquote"><blockquote class="blockquote"><p> 693 <span class="serif_italic">f<sub>1</sub>(-ζ) = ½x(-ζ) {h(-ζ)}² ⋅ b<sub>0</sub>(-ζ)</span> 694 </p></blockquote></div> 695<p> 696 where 697 </p> 698<div class="blockquote"><blockquote class="blockquote"><p> 699 <span class="serif_italic">h(-ζ) = {4(-ζ) / (x² - 1)}<sup>4</sup></span> 700 </p></blockquote></div> 701<p> 702 and 703 </p> 704<div class="blockquote"><blockquote class="blockquote"><p> 705 <span class="serif_italic">b<sub>0</sub>(-ζ) = -5/(48ζ²) + 1/(-ζ)<sup>½</sup> ⋅ { 5/(24(x<sup>2</sup>-1)<sup>3/2</sup>) + 706 1/(8(x<sup>2</sup>-1)<sup>½)</sup>}</span> 707 </p></blockquote></div> 708<p> 709 When solving for <span class="emphasis"><em>x(-ζ)</em></span> in Eq. 7 above, the right-hand-side 710 is expanded to order 2 in a Taylor series for large <span class="emphasis"><em>x</em></span>. 711 This results in 712 </p> 713<div class="blockquote"><blockquote class="blockquote"><p> 714 <span class="serif_italic">⅔(-ζ)<sup>3/2</sup> ≈ x + 1/2x - π/2</span> 715 </p></blockquote></div> 716<p> 717 The positive root of the resulting quadratic equation is used to find an 718 initial estimate <span class="emphasis"><em>x(-ζ)</em></span>. This initial estimate is subsequently 719 refined with several steps of Newton-Raphson iteration in Eq. 7. 720 </p> 721<p> 722 Estimates of the roots of cylindrical Bessel functions of negative order 723 on the positive real axis are found using interlacing relations. For example, 724 the <span class="emphasis"><em>m<sup>th</sup></em></span> root of the cylindrical Bessel function <span class="emphasis"><em>j<sub>-ν,m</sub></em></span> 725 is bracketed by the <span class="emphasis"><em>m<sup>th</sup></em></span> root and the <span class="emphasis"><em>(m+1)<sup>th</sup></em></span> 726 root of the Bessel function of corresponding positive integer order. In other 727 words, 728 </p> 729<div class="blockquote"><blockquote class="blockquote"><p> 730 <span class="serif_italic">j<sub>nν, m</sub> < j<sub>-ν, m</sub> < j<sub>nν, m+1</sub></span> 731 </p></blockquote></div> 732<p> 733 where <span class="emphasis"><em>m > 1</em></span> and <span class="emphasis"><em>n<sub>ν</sub></em></span> represents 734 the integral floor of the absolute value of <span class="emphasis"><em>|-ν|</em></span>. 735 </p> 736<p> 737 Similar bracketing relations are used to find estimates of the roots of Neumann 738 functions of negative order, whereby a discontinuity at every negative half-integer 739 order needs to be handled. 740 </p> 741<p> 742 Bracketing relations do not hold for the first root of cylindrical Bessel 743 functions and cylindrical Neumann functions with negative order. Therefore, 744 iterative algorithms combined with root-finding via bisection are used to 745 localize <span class="emphasis"><em>j<sub>-ν,1</sub></em></span> and <span class="emphasis"><em>y<sub>-ν,1</sub></em></span>. 746 </p> 747<h4> 748<a name="math_toolkit.bessel.bessel_root.h7"></a> 749 <span class="phrase"><a name="math_toolkit.bessel.bessel_root.testing"></a></span><a class="link" href="bessel_root.html#math_toolkit.bessel.bessel_root.testing">Testing</a> 750 </h4> 751<p> 752 The precision of evaluation of zeros was tested at 50 decimal digits using 753 <code class="computeroutput"><span class="identifier">cpp_dec_float_50</span></code> and found 754 identical with spot values computed by <a href="http://www.wolframalpha.com/" target="_top">Wolfram 755 Alpha</a>. 756 </p> 757</div> 758<table xmlns:rev="http://www.cs.rpi.edu/~gregod/boost/tools/doc/revision" width="100%"><tr> 759<td align="left"></td> 760<td align="right"><div class="copyright-footer">Copyright © 2006-2019 Nikhar 761 Agrawal, Anton Bikineev, Paul A. Bristow, Marco Guazzone, Christopher Kormanyos, 762 Hubert Holin, Bruno Lalande, John Maddock, Jeremy Murphy, Matthew Pulver, Johan 763 Råde, Gautam Sewani, Benjamin Sobotta, Nicholas Thompson, Thijs van den Berg, 764 Daryle Walker and Xiaogang Zhang<p> 765 Distributed under the Boost Software License, Version 1.0. (See accompanying 766 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>) 767 </p> 768</div></td> 769</tr></table> 770<hr> 771<div class="spirit-nav"> 772<a accesskey="p" href="bessel_first.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> 773</div> 774</body> 775</html> 776