1 // Copyright Nick Thompson 2017.
2 // Use, modification and distribution are subject to the
3 // Boost Software License, Version 1.0.
4 // (See accompanying file LICENSE_1_0.txt
5 // or copy at http://www.boost.org/LICENSE_1_0.txt)
6
7 #include <iostream>
8 #include <string>
9 #include <boost/math/constants/constants.hpp>
10 #include <boost/multiprecision/cpp_bin_float.hpp>
11 #include <boost/multiprecision/cpp_dec_float.hpp>
12 #include <boost/math/special_functions/legendre_stieltjes.hpp>
13
14 using boost::math::legendre_p;
15 using boost::math::legendre_p_zeros;
16 using boost::math::legendre_p_prime;
17 using boost::math::legendre_stieltjes;
18 using boost::multiprecision::cpp_bin_float_quad;
19 using boost::multiprecision::cpp_bin_float_100;
20 using boost::multiprecision::cpp_dec_float_100;
21
22 template<class Real>
gauss_kronrod_rule(size_t order)23 void gauss_kronrod_rule(size_t order)
24 {
25 std::cout << std::setprecision(std::numeric_limits<Real>::digits10);
26 std::cout << std::fixed;
27 auto gauss_nodes = boost::math::legendre_p_zeros<Real>(order);
28 auto E = legendre_stieltjes<Real>(order + 1);
29 std::vector<Real> gauss_weights(gauss_nodes.size(), std::numeric_limits<Real>::quiet_NaN());
30 std::vector<Real> gauss_kronrod_weights(gauss_nodes.size(), std::numeric_limits<Real>::quiet_NaN());
31 for (size_t i = 0; i < gauss_nodes.size(); ++i)
32 {
33 Real node = gauss_nodes[i];
34 Real lp = legendre_p_prime<Real>(order, node);
35 gauss_weights[i] = 2/( (1-node*node)*lp*lp);
36 // P_n(x) = (2n)!/(2^n (n!)^2) pi_n(x), where pi_n is the monic Legendre polynomial.
37 gauss_kronrod_weights[i] = gauss_weights[i] + static_cast<Real>(2)/(static_cast<Real>(order+1)*legendre_p_prime(order, node)*E(node));
38 }
39
40 std::cout << "static const std::vector<Real> gauss_nodes {\n";
41 for (auto const & node : gauss_nodes)
42 {
43 std::cout << " boost::lexical_cast<Real>(\"" << node << "\"),\n";
44 }
45 std::cout << "};\n\n";
46
47 std::cout << "static const std::vector<Real> gauss_weights {\n";
48 for (auto const & weight : gauss_weights)
49 {
50 std::cout << " boost::lexical_cast<Real>(\"" << weight << "\"),\n";
51 }
52 std::cout << "};\n\n";
53
54 std::cout << "static const std::vector<Real> gauss_kronrod_weights {\n";
55 for (auto const & weight : gauss_kronrod_weights)
56 {
57 std::cout << " boost::lexical_cast<Real>(\"" << weight << "\"),\n";
58 }
59 std::cout << "};\n\n";
60
61 auto kronrod_nodes = E.zeros();
62 std::vector<Real> kronrod_weights(kronrod_nodes.size());
63 for (size_t i = 0; i < kronrod_weights.size(); ++i)
64 {
65 Real node = kronrod_nodes[i];
66 kronrod_weights[i] = static_cast<Real>(2)/(static_cast<Real>(order+1)*legendre_p(order, node)*E.prime(node));
67 }
68
69 std::cout << "static const std::vector<Real> kronrod_nodes {\n";
70 for (auto node : kronrod_nodes)
71 {
72 std::cout << " boost::lexical_cast<Real>(\"" << node << "\"),\n";
73 }
74 std::cout << "};\n\n";
75
76 std::cout << "static const std::vector<Real> kronrod_weights {\n";
77 for (auto const & weight : kronrod_weights)
78 {
79 std::cout << " boost::lexical_cast<Real>(\"" << weight << "\"),\n";
80 }
81 std::cout << "};\n\n";
82
83 }
84
main()85 int main()
86 {
87 //std::cout << "Gauss-Kronrod 7-15 Rule:\n";
88 //gauss_kronrod_rule<cpp_dec_float_100>(7);
89
90 //std::cout << "\n\nGauss-Kronrod 10-21 Rule:\n";
91 //gauss_kronrod_rule<cpp_dec_float_100>(10);
92
93 std::cout << "\n\nGauss-Kronrod 15-31 Rule:\n";
94 gauss_kronrod_rule<cpp_dec_float_100>(15);
95 /*
96 std::cout << "\n\nGauss-Kronrod 20-41 Rule:\n";
97 gauss_kronrod_rule<cpp_dec_float_100>(20);
98
99 std::cout << "\n\nGauss-Kronrod 25-51 Rule:\n";
100 gauss_kronrod_rule<cpp_dec_float_100>(25);
101
102 std::cout << "\n\nGauss-Kronrod 30-61 Rule:\n";
103 gauss_kronrod_rule<cpp_dec_float_100>(30);
104
105 std::cout << "\n\nGauss-Kronrod 35-71 Rule:\n";
106 gauss_kronrod_rule<cpp_dec_float_100>(35);
107
108 std::cout << "\n\nGauss-Kronrod 40-81 Rule:\n";
109 gauss_kronrod_rule<cpp_dec_float_100>(40);*/
110 }
111