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