Skip to content

Commit 1c24566

Browse files
memmettmemmett
committed
Add pfasst-quadrature.hpp and appropriate tests.
1 parent 4fa799c commit 1c24566

File tree

2 files changed

+323
-0
lines changed

2 files changed

+323
-0
lines changed

‎src/pfasst-quadrature.hpp‎

Lines changed: 192 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,192 @@
1+
2+
#include <algorithm>
3+
#include <cmath>
4+
#include <complex>
5+
#include <vector>
6+
7+
using namespace std;
8+
9+
typedef unsigned int uint;
10+
11+
namespace pfasst {
12+
13+
template<typename ptype>
14+
class polynomial {
15+
vector<ptype> c;
16+
17+
public:
18+
19+
polynomial(uint n) : c(n) {
20+
fill(c.begin(), c.end(), 0.0);
21+
}
22+
23+
uint order() const {
24+
return c.size()-1;
25+
}
26+
27+
ptype& operator[](const unsigned int i) { return c.at(i); }
28+
29+
polynomial<ptype> differentiate() const {
30+
polynomial<ptype> p(c.size()-1);
31+
for (int j=1; j<c.size(); j++)
32+
p[j-1] = j * c[j];
33+
return p;
34+
}
35+
36+
polynomial<ptype> integrate() const {
37+
polynomial<ptype> p(c.size()+1);
38+
for (int j=0; j<c.size(); j++)
39+
p[j+1] = c[j] / (j+1);
40+
return p;
41+
}
42+
43+
template<typename xtype>
44+
xtype evaluate(const xtype x) const {
45+
int n = c.size()-1;
46+
xtype v = c[n];
47+
for (int j=n-1; j>=0; j--)
48+
v = x * v + c[j];
49+
return v;
50+
}
51+
52+
polynomial<ptype> normalize() const {
53+
polynomial<ptype> p(c.size());
54+
for (int j=0; j<c.size(); j++)
55+
p[j] = c[j] / c[c.size()-1];
56+
return p;
57+
}
58+
59+
vector<ptype> roots() const {
60+
uint n = c.size()-1;
61+
62+
// initial guess
63+
polynomial<complex<ptype>> z0(n), z1(n);
64+
for (int j=0; j<n; j++) {
65+
z0[j] = pow((0.4, 0.9), j);
66+
z1[j] = z0[j];
67+
}
68+
69+
// durand-kerner-weierstrass iterations
70+
polynomial<ptype> p = normalize();
71+
for (int k=0; k<100; k++) {
72+
complex<ptype> num, den;
73+
for (int i=0; i<n; i++) {
74+
num = p.evaluate(z0[i]);
75+
den = 1.0;
76+
for (int j=0; j<n; j++) {
77+
if (j == i) continue;
78+
den = den * (z0[i] - z0[j]);
79+
}
80+
z0[i] = z0[i] - num / den;
81+
}
82+
83+
// converged?
84+
ptype acc = 0.0;
85+
for (int j=0; j<n; j++)
86+
acc += abs(z0[j] - z1[j]);
87+
if (acc < 1e-24)
88+
break;
89+
90+
z1 = z0;
91+
}
92+
93+
vector<ptype> roots(n);
94+
for (int j=0; j<n; j++)
95+
roots[j] = abs(z0[j]) < 1e-12 ? 0.0 : real(z0[j]);
96+
97+
sort(roots.begin(), roots.end());
98+
return roots;
99+
}
100+
101+
static polynomial<ptype> legendre(const uint order)
102+
{
103+
if (order == 0) {
104+
polynomial<ptype> p(1);
105+
p[0] = 1.0;
106+
return p;
107+
}
108+
109+
if (order == 1) {
110+
polynomial<ptype> p(2);
111+
p[0] = 0.0;
112+
p[1] = 1.0;
113+
return p;
114+
}
115+
116+
polynomial<ptype> p0(order+1), p1(order+1), p2(order+1);
117+
p0[0] = 1.0; p1[1] = 1.0;
118+
119+
// (n + 1) P_{n+1} = (2n + 1) x P_{n} - n P_{n-1}
120+
for (int m=1; m<order; m++) {
121+
for (int j=1; j<order+1; j++)
122+
p2[j] = ( (2*m + 1) * p1[j-1] - m * p0[j] ) / (m + 1);
123+
p2[0] = - m * p0[0] / (m + 1);
124+
125+
for (int j=0; j<order+1; j++) {
126+
p0[j] = p1[j];
127+
p1[j] = p2[j];
128+
}
129+
}
130+
131+
return p2;
132+
}
133+
};
134+
135+
//#define pi 3.1415926535897932384626433832795028841971693993751
136+
137+
template<typename ntype>
138+
vector<ntype> compute_nodes(int nnodes, string qtype)
139+
{
140+
vector<ntype> nodes(nnodes);
141+
142+
if (qtype == "gauss-legendre") {
143+
auto roots = polynomial<ntype>::legendre(nnodes).roots();
144+
for (int j=0; j<nnodes; j++)
145+
nodes[j] = 0.5 * (1.0 + roots[j]);
146+
} else if (qtype == "gauss-lobatto") {
147+
auto roots = polynomial<ntype>::legendre(nnodes-1).differentiate().roots();
148+
for (int j=0; j<nnodes-2; j++)
149+
nodes[j+1] = 0.5 * (1.0 + roots[j]);
150+
nodes[0] = 0.0; nodes[nnodes-1] = 1.0;
151+
}
152+
153+
return nodes;
154+
}
155+
156+
157+
// void sdc_smat(sdc_mat *smat,
158+
// const int n, const int m, const int sign,
159+
// const sdc_dtype *dst, const sdc_dtype *src,
160+
// const int *flags, int ndst, int nsrc)
161+
// {
162+
// /* for (int n=0; n<(ndst-1)*nsrc; n++) */
163+
// /* smat[n] = 0.0; */
164+
165+
// sdc_dtype p[nsrc+1], p1[nsrc+1];
166+
// for (int i=0; i<nsrc; i++) {
167+
// if ((flags[i] & SDC_NODE_PROPER) == 0) continue;
168+
169+
// // construct interpolating polynomial coefficients
170+
// p[0] = 1.0; for (int j=1; j<nsrc+1; j++) p[j] = 0.0;
171+
// for (int m=0; m<nsrc; m++) {
172+
// if (((flags[m] & SDC_NODE_PROPER) == 0) || (m == i)) continue;
173+
// // p_{m+1}(x) = (x - x_j) * p_m(x)
174+
// p1[0] = 0.0;
175+
// for (int j=0; j<nsrc; j++) p1[j+1] = p[j];
176+
// for (int j=0; j<nsrc+1; j++) p1[j] -= p[j] * src[m];
177+
// for (int j=0; j<nsrc+1; j++) p[j] = p1[j];
178+
// }
179+
180+
// // evaluate integrals
181+
// sdc_dtype den = poly_eval(p, nsrc, src[i]);
182+
// poly_int(p, nsrc+1);
183+
// for (int j=1; j<ndst; j++) {
184+
// sdc_dtype s = poly_eval(p, nsrc, dst[j]) - poly_eval(p, nsrc, dst[j-1]);
185+
// sdc_mat_setvalue(smat, n+j-1, m+i, s / den, sign);
186+
// }
187+
// }
188+
// }
189+
190+
191+
192+
}

‎tests/test-quadrature.cpp‎

Lines changed: 131 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,131 @@
1+
/*
2+
* Tests for quadrature related routines: polynomials, nodes, and matrices.
3+
*/
4+
5+
#include <iostream>
6+
#include <tuple>
7+
8+
#include <gtest/gtest.h>
9+
#include <gmock/gmock.h>
10+
11+
#include <pfasst-quadrature.hpp>
12+
13+
using namespace std;
14+
15+
TEST(PolyTest, LegendrePolys) {
16+
auto l0 = pfasst::polynomial<double>::legendre(0);
17+
EXPECT_EQ(l0.order(), 0);
18+
EXPECT_EQ(l0[0], 1.0);
19+
20+
auto l1 = pfasst::polynomial<double>::legendre(1);
21+
EXPECT_EQ(l1.order(), 1);
22+
EXPECT_EQ(l1[0], 0.0);
23+
EXPECT_EQ(l1[1], 1.0);
24+
25+
auto l2 = pfasst::polynomial<double>::legendre(2);
26+
EXPECT_EQ(l2.order(), 2);
27+
EXPECT_EQ(l2[0], -0.5);
28+
EXPECT_EQ(l2[1], 0.0);
29+
EXPECT_EQ(l2[2], 1.5);
30+
31+
auto l2d = l2.differentiate();
32+
EXPECT_EQ(l2d.order(), 1);
33+
EXPECT_EQ(l2d[0], 0.0);
34+
EXPECT_EQ(l2d[1], 3.0);
35+
36+
auto l2i = l2.integrate();
37+
EXPECT_EQ(l2i.order(), 3);
38+
EXPECT_EQ(l2i[0], 0.0);
39+
EXPECT_EQ(l2i[1], -0.5);
40+
EXPECT_EQ(l2i[2], 0.0);
41+
EXPECT_EQ(l2i[3], 0.5);
42+
43+
double a1 = l2.evaluate(1.0);
44+
EXPECT_EQ(a1, 1.0);
45+
}
46+
47+
MATCHER(DoubleNear, "") {
48+
return abs(get<0>(arg) - get<1>(arg)) < 1e-15;
49+
}
50+
51+
TEST(NodesTest, GaussLegendreNodes) {
52+
const long double l3e[3] = { 0.11270166537925831,
53+
0.5,
54+
0.8872983346207417 };
55+
56+
const long double l5e[5] = { 0.046910077030668004,
57+
0.23076534494715845,
58+
0.5,
59+
0.7692346550528415,
60+
0.953089922969332 };
61+
62+
const long double l7e[7] = { 0.025446043828620736,
63+
0.12923440720030277,
64+
0.2970774243113014,
65+
0.5,
66+
0.7029225756886985,
67+
0.8707655927996972,
68+
0.9745539561713793 };
69+
70+
auto l3 = pfasst::compute_nodes<long double>(3, "gauss-legendre");
71+
EXPECT_THAT(l3, testing::Pointwise(DoubleNear(), l3e));
72+
73+
auto l5 = pfasst::compute_nodes<long double>(5, "gauss-legendre");
74+
EXPECT_THAT(l5, testing::Pointwise(DoubleNear(), l5e));
75+
76+
auto l7 = pfasst::compute_nodes<long double>(7, "gauss-legendre");
77+
EXPECT_THAT(l7, testing::Pointwise(DoubleNear(), l7e));
78+
}
79+
80+
TEST(NodesTest, GaussLobattoNodes) {
81+
const long double l2e[2] = { 0.0,
82+
1.0 };
83+
84+
const long double l3e[3] = { 0.0,
85+
0.5,
86+
1.0 };
87+
88+
const long double l5e[5] = { 0.0,
89+
0.17267316464601143,
90+
0.5,
91+
0.8273268353539885,
92+
1.0 };
93+
94+
const long double l7e[7] = { 0.0,
95+
0.08488805186071653,
96+
0.2655756032646429,
97+
0.5,
98+
0.7344243967353571,
99+
0.9151119481392834,
100+
1.0 };
101+
102+
const long double l9e[9] = { 0.0,
103+
0.05012100229426992,
104+
0.16140686024463113,
105+
0.3184412680869109,
106+
0.5,
107+
0.6815587319130891,
108+
0.8385931397553689,
109+
0.94987899770573,
110+
1.0 };
111+
112+
auto l2 = pfasst::compute_nodes<long double>(2, "gauss-lobatto");
113+
EXPECT_THAT(l2, testing::Pointwise(DoubleNear(), l2e));
114+
115+
auto l3 = pfasst::compute_nodes<long double>(3, "gauss-lobatto");
116+
EXPECT_THAT(l3, testing::Pointwise(DoubleNear(), l3e));
117+
118+
auto l5 = pfasst::compute_nodes<long double>(5, "gauss-lobatto");
119+
EXPECT_THAT(l5, testing::Pointwise(DoubleNear(), l5e));
120+
121+
auto l7 = pfasst::compute_nodes<long double>(7, "gauss-lobatto");
122+
EXPECT_THAT(l7, testing::Pointwise(DoubleNear(), l7e));
123+
124+
auto l9 = pfasst::compute_nodes<long double>(9, "gauss-lobatto");
125+
EXPECT_THAT(l9, testing::Pointwise(DoubleNear(), l9e));
126+
}
127+
128+
int main(int argc, char **argv) {
129+
testing::InitGoogleTest(&argc, argv);
130+
return RUN_ALL_TESTS();
131+
}

0 commit comments

Comments
 (0)