barycentric_interpolation_example_2.cpp 3.6 KB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108
  1. // Copyright Nick Thompson, 2017
  2. // Distributed under the Boost Software License, Version 1.0.
  3. // (See accompanying file LICENSE_1_0.txt or
  4. // copy at http://www.boost.org/LICENSE_1_0.txt).
  5. #include <iostream>
  6. #include <limits>
  7. #include <map>
  8. //[barycentric_rational_example2
  9. /*`This further example shows how to use the iterator based constructor, and then uses the
  10. function object in our root finding algorithms to locate the points where the potential
  11. achieves a specific value.
  12. */
  13. #include <boost/math/interpolators/barycentric_rational.hpp>
  14. #include <boost/range/adaptors.hpp>
  15. #include <boost/math/tools/roots.hpp>
  16. int main()
  17. {
  18. // The lithium potential is given in Kohn's paper, Table I.
  19. // (We could equally easily use an unordered_map, a list of tuples or pairs, or a 2-dimentional array).
  20. std::map<double, double> r;
  21. r[0.02] = 5.727;
  22. r[0.04] = 5.544;
  23. r[0.06] = 5.450;
  24. r[0.08] = 5.351;
  25. r[0.10] = 5.253;
  26. r[0.12] = 5.157;
  27. r[0.14] = 5.058;
  28. r[0.16] = 4.960;
  29. r[0.18] = 4.862;
  30. r[0.20] = 4.762;
  31. r[0.24] = 4.563;
  32. r[0.28] = 4.360;
  33. r[0.32] = 4.1584;
  34. r[0.36] = 3.9463;
  35. r[0.40] = 3.7360;
  36. r[0.44] = 3.5429;
  37. r[0.48] = 3.3797;
  38. r[0.52] = 3.2417;
  39. r[0.56] = 3.1209;
  40. r[0.60] = 3.0138;
  41. r[0.68] = 2.8342;
  42. r[0.76] = 2.6881;
  43. r[0.84] = 2.5662;
  44. r[0.92] = 2.4242;
  45. r[1.00] = 2.3766;
  46. r[1.08] = 2.3058;
  47. r[1.16] = 2.2458;
  48. r[1.24] = 2.2035;
  49. r[1.32] = 2.1661;
  50. r[1.40] = 2.1350;
  51. r[1.48] = 2.1090;
  52. r[1.64] = 2.0697;
  53. r[1.80] = 2.0466;
  54. r[1.96] = 2.0325;
  55. r[2.12] = 2.0288;
  56. r[2.28] = 2.0292;
  57. r[2.44] = 2.0228;
  58. r[2.60] = 2.0124;
  59. r[2.76] = 2.0065;
  60. r[2.92] = 2.0031;
  61. r[3.08] = 2.0015;
  62. r[3.24] = 2.0008;
  63. r[3.40] = 2.0004;
  64. r[3.56] = 2.0002;
  65. r[3.72] = 2.0001;
  66. // Let's discover the absissa that will generate a potential of exactly 3.0,
  67. // start by creating 2 ranges for the x and y values:
  68. auto x_range = boost::adaptors::keys(r);
  69. auto y_range = boost::adaptors::values(r);
  70. boost::math::barycentric_rational<double> b(x_range.begin(), x_range.end(), y_range.begin());
  71. //
  72. // We'll use a lamda expression to provide the functor to our root finder, since we want
  73. // the abscissa value that yields 3, not zero. We pass the functor b by value to the
  74. // lambda expression since barycentric_rational is trivial to copy.
  75. // Here we're using simple bisection to find the root:
  76. boost::uintmax_t iterations = (std::numeric_limits<boost::uintmax_t>::max)();
  77. double abscissa_3 = boost::math::tools::bisect([=](double x) { return b(x) - 3; }, 0.44, 1.24, boost::math::tools::eps_tolerance<double>(), iterations).first;
  78. std::cout << "Abscissa value that yields a potential of 3 = " << abscissa_3 << std::endl;
  79. std::cout << "Root was found in " << iterations << " iterations." << std::endl;
  80. //
  81. // However, we have a more efficient root finding algorithm than simple bisection:
  82. iterations = (std::numeric_limits<boost::uintmax_t>::max)();
  83. abscissa_3 = boost::math::tools::bracket_and_solve_root([=](double x) { return b(x) - 3; }, 0.6, 1.2, false, boost::math::tools::eps_tolerance<double>(), iterations).first;
  84. std::cout << "Abscissa value that yields a potential of 3 = " << abscissa_3 << std::endl;
  85. std::cout << "Root was found in " << iterations << " iterations." << std::endl;
  86. }
  87. //] [/barycentric_rational_example2]
  88. //[barycentric_rational_example2_out
  89. /*` Program output is:
  90. [pre
  91. Abscissa value that yields a potential of 3 = 0.604728
  92. Root was found in 54 iterations.
  93. Abscissa value that yields a potential of 3 = 0.604728
  94. Root was found in 10 iterations.
  95. ]
  96. */
  97. //]