]>
Commit | Line | Data |
---|---|---|
b32b8144 FG |
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 | #ifndef BOOST_MATH_SPECIAL_LEGENDRE_STIELTJES_HPP | |
8 | #define BOOST_MATH_SPECIAL_LEGENDRE_STIELTJES_HPP | |
9 | ||
10 | /* | |
11 | * Constructs the Legendre-Stieltjes polynomial of degree m. | |
12 | * The Legendre-Stieltjes polynomials are used to create extensions for Gaussian quadratures, | |
13 | * commonly called "Gauss-Konrod" quadratures. | |
14 | * | |
15 | * References: | |
16 | * Patterson, TNL. "The optimum addition of points to quadrature formulae." Mathematics of Computation 22.104 (1968): 847-856. | |
17 | */ | |
18 | ||
19 | #include <iostream> | |
20 | #include <vector> | |
21 | #include <boost/math/tools/roots.hpp> | |
22 | #include <boost/math/special_functions/legendre.hpp> | |
23 | ||
24 | namespace boost{ | |
25 | namespace math{ | |
26 | ||
27 | template<class Real> | |
28 | class legendre_stieltjes | |
29 | { | |
30 | public: | |
31 | legendre_stieltjes(size_t m) | |
32 | { | |
33 | if (m == 0) | |
34 | { | |
35 | throw std::domain_error("The Legendre-Stieltjes polynomial is defined for order m > 0.\n"); | |
36 | } | |
37 | m_m = static_cast<int>(m); | |
38 | std::ptrdiff_t n = m - 1; | |
39 | std::ptrdiff_t q; | |
40 | std::ptrdiff_t r; | |
41 | bool odd = n & 1; | |
42 | if (odd) | |
43 | { | |
44 | q = 1; | |
45 | r = (n-1)/2 + 2; | |
46 | } | |
47 | else | |
48 | { | |
49 | q = 0; | |
50 | r = n/2 + 1; | |
51 | } | |
52 | m_a.resize(r + 1); | |
53 | // We'll keep the ones-based indexing at the cost of storing a superfluous element | |
54 | // so that we can follow Patterson's notation exactly. | |
55 | m_a[r] = static_cast<Real>(1); | |
56 | // Make sure using the zero index is a bug: | |
57 | m_a[0] = std::numeric_limits<Real>::quiet_NaN(); | |
58 | ||
59 | for (std::ptrdiff_t k = 1; k < r; ++k) | |
60 | { | |
61 | Real ratio = 1; | |
62 | m_a[r - k] = 0; | |
63 | for (std::ptrdiff_t i = r + 1 - k; i <= r; ++i) | |
64 | { | |
65 | // See Patterson, equation 12 | |
66 | std::ptrdiff_t num = (n - q + 2*(i + k - 1))*(n + q + 2*(k - i + 1))*(n-1-q+2*(i-k))*(2*(k+i-1) -1 -q -n); | |
67 | std::ptrdiff_t den = (n - q + 2*(i - k))*(2*(k + i - 1) - q - n)*(n + 1 + q + 2*(k - i))*(n - 1 - q + 2*(i + k)); | |
68 | ratio *= static_cast<Real>(num)/static_cast<Real>(den); | |
69 | m_a[r - k] -= ratio*m_a[i]; | |
70 | } | |
71 | } | |
72 | } | |
73 | ||
74 | ||
75 | Real norm_sq() const | |
76 | { | |
77 | Real t = 0; | |
78 | bool odd = m_m & 1; | |
79 | for (size_t i = 1; i < m_a.size(); ++i) | |
80 | { | |
81 | if(odd) | |
82 | { | |
83 | t += 2*m_a[i]*m_a[i]/static_cast<Real>(4*i-1); | |
84 | } | |
85 | else | |
86 | { | |
87 | t += 2*m_a[i]*m_a[i]/static_cast<Real>(4*i-3); | |
88 | } | |
89 | } | |
90 | return t; | |
91 | } | |
92 | ||
93 | ||
94 | Real operator()(Real x) const | |
95 | { | |
96 | // Trivial implementation: | |
97 | // Em += m_a[i]*legendre_p(2*i - 1, x); m odd | |
98 | // Em += m_a[i]*legendre_p(2*i - 2, x); m even | |
99 | size_t r = m_a.size() - 1; | |
100 | Real p0 = 1; | |
101 | Real p1 = x; | |
102 | ||
103 | Real Em; | |
104 | bool odd = m_m & 1; | |
105 | if (odd) | |
106 | { | |
107 | Em = m_a[1]*p1; | |
108 | } | |
109 | else | |
110 | { | |
111 | Em = m_a[1]*p0; | |
112 | } | |
113 | ||
114 | unsigned n = 1; | |
115 | for (size_t i = 2; i <= r; ++i) | |
116 | { | |
117 | std::swap(p0, p1); | |
118 | p1 = boost::math::legendre_next(n, x, p0, p1); | |
119 | ++n; | |
120 | if (!odd) | |
121 | { | |
122 | Em += m_a[i]*p1; | |
123 | } | |
124 | std::swap(p0, p1); | |
125 | p1 = boost::math::legendre_next(n, x, p0, p1); | |
126 | ++n; | |
127 | if(odd) | |
128 | { | |
129 | Em += m_a[i]*p1; | |
130 | } | |
131 | } | |
132 | return Em; | |
133 | } | |
134 | ||
135 | ||
136 | Real prime(Real x) const | |
137 | { | |
138 | Real Em_prime = 0; | |
139 | ||
140 | for (size_t i = 1; i < m_a.size(); ++i) | |
141 | { | |
142 | if(m_m & 1) | |
143 | { | |
144 | Em_prime += m_a[i]*detail::legendre_p_prime_imp(static_cast<unsigned>(2*i - 1), x, policies::policy<>()); | |
145 | } | |
146 | else | |
147 | { | |
148 | Em_prime += m_a[i]*detail::legendre_p_prime_imp(static_cast<unsigned>(2*i - 2), x, policies::policy<>()); | |
149 | } | |
150 | } | |
151 | return Em_prime; | |
152 | } | |
153 | ||
154 | std::vector<Real> zeros() const | |
155 | { | |
156 | using boost::math::constants::half; | |
157 | ||
158 | std::vector<Real> stieltjes_zeros; | |
159 | std::vector<Real> legendre_zeros = legendre_p_zeros<Real>(m_m - 1); | |
160 | int k; | |
161 | if (m_m & 1) | |
162 | { | |
163 | stieltjes_zeros.resize(legendre_zeros.size() + 1, std::numeric_limits<Real>::quiet_NaN()); | |
164 | stieltjes_zeros[0] = 0; | |
165 | k = 1; | |
166 | } | |
167 | else | |
168 | { | |
169 | stieltjes_zeros.resize(legendre_zeros.size(), std::numeric_limits<Real>::quiet_NaN()); | |
170 | k = 0; | |
171 | } | |
172 | ||
173 | while (k < (int)stieltjes_zeros.size()) | |
174 | { | |
175 | Real lower_bound; | |
176 | Real upper_bound; | |
177 | if (m_m & 1) | |
178 | { | |
179 | lower_bound = legendre_zeros[k - 1]; | |
180 | if (k == (int)legendre_zeros.size()) | |
181 | { | |
182 | upper_bound = 1; | |
183 | } | |
184 | else | |
185 | { | |
186 | upper_bound = legendre_zeros[k]; | |
187 | } | |
188 | } | |
189 | else | |
190 | { | |
191 | lower_bound = legendre_zeros[k]; | |
192 | if (k == (int)legendre_zeros.size() - 1) | |
193 | { | |
194 | upper_bound = 1; | |
195 | } | |
196 | else | |
197 | { | |
198 | upper_bound = legendre_zeros[k+1]; | |
199 | } | |
200 | } | |
201 | ||
202 | // The root bracketing is not very tight; to keep weird stuff from happening | |
203 | // in the Newton's method, let's tighten up the tolerance using a few bisections. | |
204 | boost::math::tools::eps_tolerance<Real> tol(6); | |
205 | auto g = [&](Real t) { return this->operator()(t); }; | |
206 | auto p = boost::math::tools::bisect(g, lower_bound, upper_bound, tol); | |
207 | ||
208 | Real x_nk_guess = p.first + (p.second - p.first)*half<Real>(); | |
209 | boost::uintmax_t number_of_iterations = 500; | |
210 | ||
211 | auto f = [&] (Real x) { Real Pn = this->operator()(x); | |
212 | Real Pn_prime = this->prime(x); | |
213 | return std::pair<Real, Real>(Pn, Pn_prime); }; | |
214 | ||
215 | const Real x_nk = boost::math::tools::newton_raphson_iterate(f, x_nk_guess, | |
216 | p.first, p.second, | |
20effc67 | 217 | tools::digits<Real>(), |
b32b8144 FG |
218 | number_of_iterations); |
219 | ||
220 | BOOST_ASSERT(p.first < x_nk); | |
221 | BOOST_ASSERT(x_nk < p.second); | |
222 | stieltjes_zeros[k] = x_nk; | |
223 | ++k; | |
224 | } | |
225 | return stieltjes_zeros; | |
226 | } | |
227 | ||
228 | private: | |
229 | // Coefficients of Legendre expansion | |
230 | std::vector<Real> m_a; | |
231 | int m_m; | |
232 | }; | |
233 | ||
234 | }} | |
235 | #endif |