]>
Commit | Line | Data |
---|---|---|
92f5a8d4 TL |
1 | /* |
2 | * Copyright Nick Thompson, 2019 | |
3 | * Use, modification and distribution are subject to the | |
4 | * Boost Software License, Version 1.0. (See accompanying file | |
5 | * LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt) | |
6 | */ | |
7 | ||
8 | #include "math_unit_test.hpp" | |
9 | #include <numeric> | |
10 | #include <utility> | |
11 | #include <boost/math/interpolators/cardinal_quadratic_b_spline.hpp> | |
12 | using boost::math::interpolators::cardinal_quadratic_b_spline; | |
13 | ||
14 | template<class Real> | |
15 | void test_constant() | |
16 | { | |
17 | Real c = 7.2; | |
18 | Real t0 = 0; | |
19 | Real h = Real(1)/Real(16); | |
20 | size_t n = 512; | |
21 | std::vector<Real> v(n, c); | |
22 | auto qbs = cardinal_quadratic_b_spline<Real>(v.data(), v.size(), t0, h); | |
23 | ||
24 | size_t i = 0; | |
25 | while (i < n) { | |
26 | Real t = t0 + i*h; | |
27 | CHECK_ULP_CLOSE(c, qbs(t), 2); | |
28 | CHECK_MOLLIFIED_CLOSE(0, qbs.prime(t), 100*std::numeric_limits<Real>::epsilon()); | |
29 | ++i; | |
30 | } | |
31 | ||
32 | i = 0; | |
33 | while (i < n) { | |
34 | Real t = t0 + i*h + h/2; | |
35 | CHECK_ULP_CLOSE(c, qbs(t), 2); | |
36 | CHECK_MOLLIFIED_CLOSE(0, qbs.prime(t), 300*std::numeric_limits<Real>::epsilon()); | |
37 | t = t0 + i*h + h/4; | |
38 | CHECK_ULP_CLOSE(c, qbs(t), 2); | |
39 | CHECK_MOLLIFIED_CLOSE(0, qbs.prime(t), 150*std::numeric_limits<Real>::epsilon()); | |
40 | ++i; | |
41 | } | |
42 | } | |
43 | ||
44 | template<class Real> | |
45 | void test_linear() | |
46 | { | |
47 | Real m = 8.3; | |
48 | Real b = 7.2; | |
49 | Real t0 = 0; | |
50 | Real h = Real(1)/Real(16); | |
51 | size_t n = 512; | |
52 | std::vector<Real> y(n); | |
53 | for (size_t i = 0; i < n; ++i) { | |
54 | Real t = i*h; | |
55 | y[i] = m*t + b; | |
56 | } | |
57 | auto qbs = cardinal_quadratic_b_spline<Real>(y.data(), y.size(), t0, h); | |
58 | ||
59 | size_t i = 0; | |
60 | while (i < n) { | |
61 | Real t = t0 + i*h; | |
62 | CHECK_ULP_CLOSE(m*t+b, qbs(t), 2); | |
63 | CHECK_ULP_CLOSE(m, qbs.prime(t), 820); | |
64 | ++i; | |
65 | } | |
66 | ||
67 | i = 0; | |
68 | while (i < n) { | |
69 | Real t = t0 + i*h + h/2; | |
70 | CHECK_ULP_CLOSE(m*t+b, qbs(t), 2); | |
71 | CHECK_MOLLIFIED_CLOSE(m, qbs.prime(t), 1500*std::numeric_limits<Real>::epsilon()); | |
72 | t = t0 + i*h + h/4; | |
73 | CHECK_ULP_CLOSE(m*t+b, qbs(t), 3); | |
74 | CHECK_MOLLIFIED_CLOSE(m, qbs.prime(t), 1500*std::numeric_limits<Real>::epsilon()); | |
75 | ++i; | |
76 | } | |
77 | } | |
78 | ||
79 | template<class Real> | |
80 | void test_quadratic() | |
81 | { | |
82 | Real a = 8.2; | |
83 | Real b = 7.2; | |
84 | Real c = -9.2; | |
85 | Real t0 = 0; | |
86 | Real h = Real(1)/Real(16); | |
87 | size_t n = 513; | |
88 | std::vector<Real> y(n); | |
89 | for (size_t i = 0; i < n; ++i) { | |
90 | Real t = i*h; | |
91 | y[i] = a*t*t + b*t + c; | |
92 | } | |
93 | Real t_max = t0 + (n-1)*h; | |
94 | auto qbs = cardinal_quadratic_b_spline<Real>(y, t0, h, b, 2*a*t_max + b); | |
95 | ||
96 | size_t i = 0; | |
97 | while (i < n) { | |
98 | Real t = t0 + i*h; | |
99 | CHECK_ULP_CLOSE(a*t*t + b*t + c, qbs(t), 2); | |
100 | ++i; | |
101 | } | |
102 | ||
103 | i = 0; | |
104 | while (i < n) { | |
105 | Real t = t0 + i*h + h/2; | |
106 | CHECK_ULP_CLOSE(a*t*t + b*t + c, qbs(t), 47); | |
107 | ||
108 | t = t0 + i*h + h/4; | |
109 | if (!CHECK_ULP_CLOSE(a*t*t + b*t + c, qbs(t), 104)) { | |
110 | std::cerr << " Problem abscissa t = " << t << "\n"; | |
111 | } | |
112 | ++i; | |
113 | } | |
114 | } | |
115 | ||
116 | int main() | |
117 | { | |
118 | test_constant<float>(); | |
119 | test_constant<double>(); | |
1e59de90 | 120 | #ifndef BOOST_MATH_NO_LONG_DOUBLE_MATH_FUNCTIONS |
92f5a8d4 | 121 | test_constant<long double>(); |
1e59de90 | 122 | #endif |
92f5a8d4 TL |
123 | |
124 | test_linear<float>(); | |
125 | test_linear<double>(); | |
1e59de90 | 126 | #ifndef BOOST_MATH_NO_LONG_DOUBLE_MATH_FUNCTIONS |
92f5a8d4 | 127 | test_linear<long double>(); |
1e59de90 | 128 | #endif |
92f5a8d4 TL |
129 | |
130 | test_quadratic<double>(); | |
1e59de90 | 131 | #ifndef BOOST_MATH_NO_LONG_DOUBLE_MATH_FUNCTIONS |
92f5a8d4 | 132 | test_quadratic<long double>(); |
1e59de90 | 133 | #endif |
92f5a8d4 TL |
134 | return boost::math::test::report_errors(); |
135 | } |