ceph/src/boost/libs/math/example/bessel_zeros_example.cpp

   1 // Copyright Christopher Kormanyos 2013.
   2 // Copyright Paul A. Bristow 2013.
   3 // Copyright John Maddock 2013.
   4
   5 // Distributed under the Boost Software License, Version 1.0.
   6 // (See accompanying file LICENSE_1_0.txt or
   7 // copy at http://www.boost.org/LICENSE_1_0.txt).
   8
   9 #ifdef _MSC_VER
  10 #  pragma warning (disable : 4512) // assignment operator could not be generated.
  11 #  pragma warning (disable : 4996) // assignment operator could not be generated.
  12 #endif
  13
  14 #include <iostream>
  15 #include <limits>
  16 #include <vector>
  17 #include <algorithm>
  18 #include <iomanip>
  19 #include <iterator>
  20
  21 // Weisstein, Eric W. "Bessel Function Zeros." From MathWorld--A Wolfram Web Resource.
  22 // http://mathworld.wolfram.com/BesselFunctionZeros.html
  23 // Test values can be calculated using [@wolframalpha.com WolframAplha]
  24 // See also http://dlmf.nist.gov/10.21
  25
  26 //[bessel_zero_example_1
  27
  28 /*`This example demonstrates calculating zeros of the Bessel, Neumann and Airy functions.
  29 It also shows how Boost.Math and Boost.Multiprecision can be combined to provide
  30 a many decimal digit precision. For 50 decimal digit precision we need to include
  31 */
  32
  33   #include <boost/multiprecision/cpp_dec_float.hpp>
  34
  35 /*`and a `typedef` for `float_type` may be convenient
  36 (allowing a quick switch to re-compute at built-in `double` or other precision)
  37 */
  38   typedef boost::multiprecision::cpp_dec_float_50 float_type;
  39
  40 //`To use the functions for finding zeros of the functions we need
  41
  42   #include <boost/math/special_functions/bessel.hpp>
  43
  44 //`This file includes the forward declaration signatures for the zero-finding functions:
  45
  46 //  #include <boost/math/special_functions/math_fwd.hpp>
  47
  48 /*`but more details are in the full documentation, for example at
  49 [@http://www.boost.org/doc/libs/1_53_0/libs/math/doc/sf_and_dist/html/math_toolkit/special/bessel/bessel_over.html Boost.Math Bessel functions]
  50 */
  51
  52 /*`This example shows obtaining both a single zero of the Bessel function,
  53 and then placing multiple zeros into a container like `std::vector` by providing an iterator.
  54 The signature of the single value function is:
  55
  56   template <class T>
  57   inline typename detail::bessel_traits<T, T, policies::policy<> >::result_type
  58     cyl_bessel_j_zero(T v,  // Floating-point value for Jv.
  59     int m); // start index.
  60
  61 The result type is controlled by the floating-point type of parameter `v`
  62 (but subject to the usual __precision_policy and __promotion_policy).
  63
  64 The signature of multiple zeros function is:
  65
  66   template <class T, class OutputIterator>
  67   inline OutputIterator cyl_bessel_j_zero(T v, // Floating-point value for Jv.
  68                                 int start_index, // 1-based start index.
  69                                 unsigned number_of_zeros,
  70                                 OutputIterator out_it); // iterator into container for zeros.
  71
  72 There is also a version which allows control of the __policy_section for error handling and precision.
  73
  74   template <class T, class OutputIterator, class Policy>
  75   inline OutputIterator cyl_bessel_j_zero(T v, // Floating-point value for Jv.
  76                                 int start_index, // 1-based start index.
  77                                 unsigned number_of_zeros,
  78                                 OutputIterator out_it,
  79                                 const Policy& pol); // iterator into container for zeros.
  80
  81 */
  82 //]  [/bessel_zero_example_1]
  83
  84 //[bessel_zero_example_iterator_1]
  85 /*`We use the `cyl_bessel_j_zero` output iterator parameter `out_it`
  86 to create a sum of 1/zeros[super 2] by defining a custom output iterator:
  87 */
  88
  89 template <class T>
  90 struct output_summation_iterator
  91 {
  92    output_summation_iterator(T* p) : p_sum(p)
  93    {}
  94    output_summation_iterator& operator*()
  95    { return *this; }
  96     output_summation_iterator& operator++()
  97    { return *this; }
  98    output_summation_iterator& operator++(int)
  99    { return *this; }
 100    output_summation_iterator& operator = (T const& val)
 101    {
 102      *p_sum += 1./ (val * val); // Summing 1/zero^2.
 103      return *this;
 104    }
 105 private:
 106    T* p_sum;
 107 };
 108
 109
 110 //] [/bessel_zero_example_iterator_1]
 111
 112 int main()
 113 {
 114   try
 115   {
 116 //[bessel_zero_example_2]
 117
 118 /*`[tip It is always wise to place code using Boost.Math inside try'n'catch blocks;
 119 this will ensure that helpful error messages can be shown when exceptional conditions arise.]
 120
 121 First, evaluate a single Bessel zero.
 122
 123 The precision is controlled by the float-point type of template parameter `T` of `v`
 124 so this example has `double` precision, at least 15 but up to 17 decimal digits (for the common 64-bit double).
 125 */
 126     double root = boost::math::cyl_bessel_j_zero(0.0, 1);
 127     // Displaying with default precision of 6 decimal digits:
 128     std::cout << "boost::math::cyl_bessel_j_zero(0.0, 1) " << root << std::endl; // 2.40483
 129     // And with all the guaranteed (15) digits:
 130     std::cout.precision(std::numeric_limits<double>::digits10);
 131     std::cout << "boost::math::cyl_bessel_j_zero(0.0, 1) " << root << std::endl; // 2.40482555769577
 132 /*`But note that because the parameter `v` controls the precision of the result,
 133 `v` [*must be a floating-point type].
 134 So if you provide an integer type, say 0, rather than 0.0, then it will fail to compile thus:
 135 ``
 136     root = boost::math::cyl_bessel_j_zero(0, 1);
 137 ``
 138 with this error message
 139 ``
 140   error C2338: Order must be a floating-point type.
 141 ``
 142
 143 Optionally, we can use a policy to ignore errors, C-style, returning some value
 144 perhaps infinity or NaN, or the best that can be done. (See __user_error_handling).
 145
 146 To create a (possibly unwise!) policy that ignores all errors:
 147 */
 148
 149   typedef boost::math::policies::policy
 150     <
 151       boost::math::policies::domain_error<boost::math::policies::ignore_error>,
 152       boost::math::policies::overflow_error<boost::math::policies::ignore_error>,
 153       boost::math::policies::underflow_error<boost::math::policies::ignore_error>,
 154       boost::math::policies::denorm_error<boost::math::policies::ignore_error>,
 155       boost::math::policies::pole_error<boost::math::policies::ignore_error>,
 156       boost::math::policies::evaluation_error<boost::math::policies::ignore_error>
 157     > ignore_all_policy;
 158
 159     double inf = std::numeric_limits<double>::infinity();
 160     double nan = std::numeric_limits<double>::quiet_NaN();
 161
 162     std::cout << "boost::math::cyl_bessel_j_zero(-1.0, 0) " << std::endl;
 163     double dodgy_root = boost::math::cyl_bessel_j_zero(-1.0, 0, ignore_all_policy());
 164     std::cout << "boost::math::cyl_bessel_j_zero(-1.0, 1) " << dodgy_root << std::endl; // 1.#QNAN
 165     double inf_root = boost::math::cyl_bessel_j_zero(inf, 1, ignore_all_policy());
 166     std::cout << "boost::math::cyl_bessel_j_zero(inf, 1) " << inf_root << std::endl; // 1.#QNAN
 167     double nan_root = boost::math::cyl_bessel_j_zero(nan, 1, ignore_all_policy());
 168     std::cout << "boost::math::cyl_bessel_j_zero(nan, 1) " << nan_root << std::endl; // 1.#QNAN
 169
 170 /*`Another version of `cyl_bessel_j_zero` allows calculation of multiple zeros with one call,
 171 placing the results in a container, often `std::vector`.
 172 For example, generate five `double` roots of J[sub v] for integral order 2.
 173
 174 showing the same results as column J[sub 2](x) in table 1 of
 175 [@ http://mathworld.wolfram.com/BesselFunctionZeros.html Wolfram Bessel Function Zeros].
 176
 177 */
 178     unsigned int n_roots = 5U;
 179     std::vector<double> roots;
 180     boost::math::cyl_bessel_j_zero(2.0, 1, n_roots, std::back_inserter(roots));
 181     std::copy(roots.begin(),
 182               roots.end(),
 183               std::ostream_iterator<double>(std::cout, "\n"));
 184
 185 /*`Or generate 50 decimal digit roots of J[sub v] for non-integral order `v = 71/19`.
 186
 187 We set the precision of the output stream and show trailing zeros to display a fixed 50 decimal digits.
 188 */
 189     std::cout.precision(std::numeric_limits<float_type>::digits10); // 50 decimal digits.
 190     std::cout << std::showpoint << std::endl; // Show trailing zeros.
 191
 192     float_type x = float_type(71) / 19;
 193     float_type r = boost::math::cyl_bessel_j_zero(x, 1); // 1st root.
 194     std::cout << "x = " << x << ", r = " << r << std::endl;
 195
 196     r = boost::math::cyl_bessel_j_zero(x, 20U); // 20th root.
 197     std::cout << "x = " << x << ", r = " << r << std::endl;
 198
 199     std::vector<float_type> zeros;
 200     boost::math::cyl_bessel_j_zero(x, 1, 3, std::back_inserter(zeros));
 201
 202     std::cout << "cyl_bessel_j_zeros" << std::endl;
 203     // Print the roots to the output stream.
 204     std::copy(zeros.begin(), zeros.end(),
 205               std::ostream_iterator<float_type>(std::cout, "\n"));
 206
 207 /*`The Neumann function zeros are evaluated very similarly:
 208 */
 209     using boost::math::cyl_neumann_zero;
 210
 211     double zn = cyl_neumann_zero(2., 1);
 212
 213     std::cout << "cyl_neumann_zero(2., 1) = " << std::endl;
 214     //double zn0 = zn;
 215     //    std::cout << "zn0 = " << std::endl;
 216     //    std::cout << zn0 << std::endl;
 217     //
 218     std::cout << zn << std::endl;
 219     //  std::cout << cyl_neumann_zero(2., 1) << std::endl;
 220
 221     std::vector<float> nzeros(3); // Space for 3 zeros.
 222     cyl_neumann_zero<float>(2.F, 1, nzeros.size(), nzeros.begin());
 223
 224     std::cout << "cyl_neumann_zero<float>(2.F, 1, " << std::endl;
 225     // Print the zeros to the output stream.
 226     std::copy(nzeros.begin(), nzeros.end(),
 227               std::ostream_iterator<float>(std::cout, "\n"));
 228
 229     std::cout << cyl_neumann_zero(static_cast<float_type>(220)/100, 1) << std::endl;
 230     // 3.6154383428745996706772556069431792744372398748422
 231
 232 /*`Finally we show how the output iterator can be used to compute a sum of zeros.
 233
 234 (See [@https://doi.org/10.1017/S2040618500034067 Ian N. Sneddon, Infinite Sums of Bessel Zeros],
 235 page 150 equation 40).
 236 */
 237 //] [/bessel_zero_example_2]
 238
 239     {
 240 //[bessel_zero_example_iterator_2]
 241 /*`The sum is calculated for many values, converging on the analytical exact value of `1/8`.
 242 */
 243     using boost::math::cyl_bessel_j_zero;
 244     double nu = 1.;
 245     double sum = 0;
 246     output_summation_iterator<double> it(&sum);  // sum of 1/zeros^2
 247     cyl_bessel_j_zero(nu, 1, 10000, it);
 248
 249     double s = 1/(4 * (nu + 1)); // 0.125 = 1/8 is exact analytical solution.
 250     std::cout << std::setprecision(6) << "nu = " << nu << ", sum = " << sum
 251       << ", exact = " << s << std::endl;
 252     // nu = 1.00000, sum = 0.124990, exact = 0.125000
 253 //] [/bessel_zero_example_iterator_2]
 254     }
 255   }
 256   catch (std::exception& ex)
 257   {
 258     std::cout << "Thrown exception " << ex.what() << std::endl;
 259   }
 260
 261 //[bessel_zero_example_iterator_3]
 262
 263 /*`Examples below show effect of 'bad' arguments that throw a `domain_error` exception.
 264 */
 265   try
 266   { // Try a negative rank m.
 267     std::cout << "boost::math::cyl_bessel_j_zero(-1.F, -1) " << std::endl;
 268     float dodgy_root = boost::math::cyl_bessel_j_zero(-1.F, -1);
 269     std::cout << "boost::math::cyl_bessel_j_zero(-1.F, -1) " << dodgy_root << std::endl;
 270     // Throw exception Error in function boost::math::cyl_bessel_j_zero<double>(double, int):
 271     // Order argument is -1, but must be >= 0 !
 272   }
 273   catch (std::exception& ex)
 274   {
 275     std::cout << "Throw exception " << ex.what() << std::endl;
 276   }
 277
 278 /*`[note The type shown is the type [*after promotion],
 279 using __precision_policy and __promotion_policy, from `float` to `double` in this case.]
 280
 281 In this example the promotion goes:
 282
 283 # Arguments are `float` and `int`.
 284 # Treat `int` "as if" it were a `double`, so arguments are `float` and `double`.
 285 # Common type is `double` - so that's the precision we want (and the type that will be returned).
 286 # Evaluate internally as `long double` for full `double` precision.
 287
 288 See full code for other examples that promote from `double` to `long double`.
 289
 290 */
 291
 292 //] [/bessel_zero_example_iterator_3]
 293     try
 294   { // order v = inf
 295      std::cout << "boost::math::cyl_bessel_j_zero(infF, 1) " << std::endl;
 296      float infF = std::numeric_limits<float>::infinity();
 297      float inf_root = boost::math::cyl_bessel_j_zero(infF, 1);
 298       std::cout << "boost::math::cyl_bessel_j_zero(infF, 1) " << inf_root << std::endl;
 299      //  boost::math::cyl_bessel_j_zero(-1.F, -1)
 300      //Thrown exception Error in function boost::math::cyl_bessel_j_zero<double>(double, int):
 301      // Requested the -1'th zero, but the rank must be positive !
 302   }
 303   catch (std::exception& ex)
 304   {
 305     std::cout << "Thrown exception " << ex.what() << std::endl;
 306   }
 307   try
 308   { // order v = inf
 309      double inf = std::numeric_limits<double>::infinity();
 310      double inf_root = boost::math::cyl_bessel_j_zero(inf, 1);
 311      std::cout << "boost::math::cyl_bessel_j_zero(inf, 1) " << inf_root << std::endl;
 312      // Throw exception Error in function boost::math::cyl_bessel_j_zero<long double>(long double, unsigned):
 313      // Order argument is 1.#INF, but must be finite >= 0 !
 314   }
 315   catch (std::exception& ex)
 316   {
 317     std::cout << "Thrown exception " << ex.what() << std::endl;
 318   }
 319
 320   try
 321   { // order v = NaN
 322      double nan = std::numeric_limits<double>::quiet_NaN();
 323      double nan_root = boost::math::cyl_bessel_j_zero(nan, 1);
 324      std::cout << "boost::math::cyl_bessel_j_zero(nan, 1) " << nan_root << std::endl;
 325      // Throw exception Error in function boost::math::cyl_bessel_j_zero<long double>(long double, unsigned):
 326      // Order argument is 1.#QNAN, but must be finite >= 0 !
 327   }
 328   catch (std::exception& ex)
 329   {
 330     std::cout << "Thrown exception " << ex.what() << std::endl;
 331   }
 332
 333   try
 334   {   // Try a negative m.
 335     double dodgy_root = boost::math::cyl_bessel_j_zero(0.0, -1);
 336     //  warning C4146: unary minus operator applied to unsigned type, result still unsigned.
 337     std::cout << "boost::math::cyl_bessel_j_zero(0.0, -1) " << dodgy_root << std::endl;
 338     //  boost::math::cyl_bessel_j_zero(0.0, -1) 6.74652e+009
 339     // This *should* fail because m is unreasonably large.
 340
 341   }
 342   catch (std::exception& ex)
 343   {
 344     std::cout << "Thrown exception " << ex.what() << std::endl;
 345   }
 346
 347   try
 348   { // m = inf
 349      double inf = std::numeric_limits<double>::infinity();
 350      double inf_root = boost::math::cyl_bessel_j_zero(0.0, inf);
 351      // warning C4244: 'argument' : conversion from 'double' to 'int', possible loss of data.
 352      std::cout << "boost::math::cyl_bessel_j_zero(0.0, inf) " << inf_root << std::endl;
 353      // Throw exception Error in function boost::math::cyl_bessel_j_zero<long double>(long double, int):
 354      // Requested the 0'th zero, but must be > 0 !
 355
 356   }
 357   catch (std::exception& ex)
 358   {
 359     std::cout << "Thrown exception " << ex.what() << std::endl;
 360   }
 361
 362   try
 363   { // m = NaN
 364      std::cout << "boost::math::cyl_bessel_j_zero(0.0, nan) " << std::endl ;
 365      double nan = std::numeric_limits<double>::quiet_NaN();
 366      double nan_root = boost::math::cyl_bessel_j_zero(0.0, nan);
 367      // warning C4244: 'argument' : conversion from 'double' to 'int', possible loss of data.
 368      std::cout << nan_root << std::endl;
 369      // Throw exception Error in function boost::math::cyl_bessel_j_zero<long double>(long double, int):
 370      // Requested the 0'th zero, but must be > 0 !
 371   }
 372   catch (std::exception& ex)
 373   {
 374     std::cout << "Thrown exception " << ex.what() << std::endl;
 375   }
 376
 377  } // int main()
 378
 379 /*
 380 Mathematica: Table[N[BesselJZero[71/19, n], 50], {n, 1, 20, 1}]
 381
 382 7.2731751938316489503185694262290765588963196701623
 383 10.724858308883141732536172745851416647110749599085
 384 14.018504599452388106120459558042660282427471931581
 385 17.25249845917041718216248716654977734919590383861
 386 20.456678874044517595180234083894285885460502077814
 387 23.64363089714234522494551422714731959985405172504
 388 26.819671140255087745421311470965019261522390519297
 389 29.988343117423674742679141796661432043878868194142
 390 33.151796897690520871250862469973445265444791966114
 391 36.3114160002162074157243540350393860813165201842
 392 39.468132467505236587945197808083337887765967032029
 393 42.622597801391236474855034831297954018844433480227
 394 45.775281464536847753390206207806726581495950012439
 395 48.926530489173566198367766817478553992471739894799
 396 52.076607045343002794279746041878924876873478063472
 397 55.225712944912571393594224327817265689059002890192
 398 58.374006101538886436775188150439025201735151418932
 399 61.521611873000965273726742659353136266390944103571
 400 64.66863105379093036834648221487366079456596628716
 401 67.815145619696290925556791375555951165111460585458
 402
 403 Mathematica: Table[N[BesselKZero[2, n], 50], {n, 1, 5, 1}]
 404 n |
 405 1 | 3.3842417671495934727014260185379031127323883259329
 406 2 | 6.7938075132682675382911671098369487124493222183854
 407 3 | 10.023477979360037978505391792081418280789658279097
 408
 409
 410 */
 411
 412  /*
 413 [bessel_zero_output]
 414
 415   boost::math::cyl_bessel_j_zero(0.0, 1) 2.40483
 416   boost::math::cyl_bessel_j_zero(0.0, 1) 2.40482555769577
 417   boost::math::cyl_bessel_j_zero(-1.0, 1) 1.#QNAN
 418   boost::math::cyl_bessel_j_zero(inf, 1) 1.#QNAN
 419   boost::math::cyl_bessel_j_zero(nan, 1) 1.#QNAN
 420   5.13562230184068
 421   8.41724414039986
 422   11.6198411721491
 423   14.7959517823513
 424   17.9598194949878
 425
 426   x = 3.7368421052631578947368421052631578947368421052632, r = 7.2731751938316489503185694262290765588963196701623
 427   x = 3.7368421052631578947368421052631578947368421052632, r = 67.815145619696290925556791375555951165111460585458
 428   7.2731751938316489503185694262290765588963196701623
 429   10.724858308883141732536172745851416647110749599085
 430   14.018504599452388106120459558042660282427471931581
 431   cyl_neumann_zero(2., 1) = 3.3842417671495935000000000000000000000000000000000
 432   3.3842418193817139000000000000000000000000000000000
 433   6.7938075065612793000000000000000000000000000000000
 434   10.023477554321289000000000000000000000000000000000
 435   3.6154383428745996706772556069431792744372398748422
 436   nu = 1.00000, sum = 0.124990, exact = 0.125000
 437   Throw exception Error in function boost::math::cyl_bessel_j_zero<double>(double, int): Order argument is -1, but must be >= 0 !
 438   Throw exception Error in function boost::math::cyl_bessel_j_zero<long double>(long double, int): Order argument is 1.#INF, but must be finite >= 0 !
 439   Throw exception Error in function boost::math::cyl_bessel_j_zero<long double>(long double, int): Order argument is 1.#QNAN, but must be finite >= 0 !
 440   Throw exception Error in function boost::math::cyl_bessel_j_zero<long double>(long double, int): Requested the -1'th zero, but must be > 0 !
 441   Throw exception Error in function boost::math::cyl_bessel_j_zero<long double>(long double, int): Requested the -2147483648'th zero, but must be > 0 !
 442   Throw exception Error in function boost::math::cyl_bessel_j_zero<long double>(long double, int): Requested the -2147483648'th zero, but must be > 0 !
 443
 444
 445 ] [/bessel_zero_output]
 446 */
 447