]>
Commit | Line | Data |
---|---|---|
b32b8144 FG |
1 | /* |
2 | * Copyright Nick Thompson, 2017 | |
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 | */ | |
20effc67 | 7 | #include "math_unit_test.hpp" |
b32b8144 | 8 | #include <boost/type_index.hpp> |
b32b8144 FG |
9 | #include <boost/math/special_functions/chebyshev.hpp> |
10 | #include <boost/math/special_functions/chebyshev_transform.hpp> | |
11 | #include <boost/math/special_functions/sinc.hpp> | |
b32b8144 FG |
12 | |
13 | #if !defined(TEST1) && !defined(TEST2) && !defined(TEST3) && !defined(TEST4) | |
14 | # define TEST1 | |
15 | # define TEST2 | |
16 | # define TEST3 | |
17 | # define TEST4 | |
18 | #endif | |
19 | ||
b32b8144 FG |
20 | using boost::math::chebyshev_t; |
21 | using boost::math::chebyshev_t_prime; | |
22 | using boost::math::chebyshev_u; | |
23 | using boost::math::chebyshev_transform; | |
24 | ||
25 | ||
26 | template<class Real> | |
27 | void test_sin_chebyshev_transform() | |
28 | { | |
29 | using boost::math::chebyshev_transform; | |
30 | using boost::math::constants::half_pi; | |
31 | using std::sin; | |
32 | using std::cos; | |
33 | using std::abs; | |
34 | ||
20effc67 TL |
35 | Real tol = std::numeric_limits<Real>::epsilon(); |
36 | auto f = [](Real x)->Real { return sin(x); }; | |
b32b8144 FG |
37 | Real a = 0; |
38 | Real b = 1; | |
39 | chebyshev_transform<Real> cheb(f, a, b, tol); | |
40 | ||
41 | Real x = a; | |
42 | while (x < b) | |
43 | { | |
44 | Real s = sin(x); | |
45 | Real c = cos(x); | |
20effc67 TL |
46 | CHECK_ABSOLUTE_ERROR(s, cheb(x), tol); |
47 | CHECK_ABSOLUTE_ERROR(c, cheb.prime(x), 150*tol); | |
b32b8144 FG |
48 | x += static_cast<Real>(1)/static_cast<Real>(1 << 7); |
49 | } | |
50 | ||
51 | Real Q = cheb.integrate(); | |
52 | ||
20effc67 | 53 | CHECK_ABSOLUTE_ERROR(1 - cos(static_cast<Real>(1)), Q, 100*tol); |
b32b8144 FG |
54 | } |
55 | ||
56 | ||
57 | template<class Real> | |
58 | void test_sinc_chebyshev_transform() | |
59 | { | |
60 | using std::cos; | |
61 | using std::sin; | |
62 | using std::abs; | |
63 | using boost::math::sinc_pi; | |
64 | using boost::math::chebyshev_transform; | |
65 | using boost::math::constants::half_pi; | |
66 | ||
20effc67 | 67 | Real tol = 100*std::numeric_limits<Real>::epsilon(); |
b32b8144 FG |
68 | auto f = [](Real x) { return boost::math::sinc_pi(x); }; |
69 | Real a = 0; | |
70 | Real b = 1; | |
71 | chebyshev_transform<Real> cheb(f, a, b, tol/50); | |
72 | ||
73 | Real x = a; | |
74 | while (x < b) | |
75 | { | |
76 | Real s = sinc_pi(x); | |
77 | Real ds = (cos(x)-sinc_pi(x))/x; | |
78 | if (x == 0) { ds = 0; } | |
b32b8144 | 79 | |
20effc67 TL |
80 | CHECK_ABSOLUTE_ERROR(s, cheb(x), tol); |
81 | CHECK_ABSOLUTE_ERROR(ds, cheb.prime(x), 10*tol); | |
b32b8144 FG |
82 | x += static_cast<Real>(1)/static_cast<Real>(1 << 7); |
83 | } | |
84 | ||
85 | Real Q = cheb.integrate(); | |
86 | //NIntegrate[Sinc[x], {x, 0, 1}, WorkingPrecision -> 200, AccuracyGoal -> 150, PrecisionGoal -> 150, MaxRecursion -> 150] | |
87 | Real Q_exp = boost::lexical_cast<Real>("0.94608307036718301494135331382317965781233795473811179047145477356668"); | |
20effc67 | 88 | CHECK_ABSOLUTE_ERROR(Q_exp, Q, tol); |
b32b8144 FG |
89 | } |
90 | ||
91 | ||
92 | ||
93 | //Examples taken from "Approximation Theory and Approximation Practice", by Trefethen | |
94 | template<class Real> | |
95 | void test_atap_examples() | |
96 | { | |
97 | using std::sin; | |
20effc67 TL |
98 | using std::exp; |
99 | using std::sqrt; | |
b32b8144 FG |
100 | using boost::math::constants::half; |
101 | using boost::math::sinc_pi; | |
102 | using boost::math::chebyshev_transform; | |
103 | using boost::math::constants::half_pi; | |
104 | ||
105 | Real tol = 10*std::numeric_limits<Real>::epsilon(); | |
106 | auto f1 = [](Real x) { return ((0 < x) - (x < 0)) - x/2; }; | |
107 | auto f2 = [](Real x) { Real t = sin(6*x); Real s = sin(x + exp(2*x)); | |
108 | Real u = (0 < s) - (s < 0); | |
109 | return t + u; }; | |
110 | ||
20effc67 TL |
111 | //auto f3 = [](Real x) { return sin(6*x) + sin(60*exp(x)); }; |
112 | //auto f4 = [](Real x) { return 1/(1+1000*(x+half<Real>())*(x+half<Real>())) + 1/sqrt(1+1000*(x-Real(1)/Real(2))*(x-Real(1)/Real(2)));}; | |
b32b8144 FG |
113 | Real a = -1; |
114 | Real b = 1; | |
20effc67 | 115 | chebyshev_transform<Real> cheb1(f1, a, b, tol); |
b32b8144 FG |
116 | chebyshev_transform<Real> cheb2(f2, a, b, tol); |
117 | //chebyshev_transform<Real> cheb3(f3, a, b, tol); | |
118 | ||
119 | Real x = a; | |
120 | while (x < b) | |
121 | { | |
20effc67 TL |
122 | // f1 and f2 are not differentiable; standard convergence rate theorems don't apply. |
123 | // Basically, the max refinements are always hit; so the error is not related to the precision of the type. | |
124 | Real acceptable_error = sqrt(tol); | |
125 | Real acceptable_error_2 = 9e-4; | |
126 | if (std::is_same<Real, long double>::value) | |
b32b8144 | 127 | { |
20effc67 | 128 | acceptable_error = 1.6e-5; |
b32b8144 | 129 | } |
20effc67 | 130 | if (std::is_same<Real, double>::value) |
b32b8144 | 131 | { |
20effc67 | 132 | acceptable_error *= 500; |
b32b8144 | 133 | } |
20effc67 TL |
134 | CHECK_ABSOLUTE_ERROR(f1(x), cheb1(x), acceptable_error); |
135 | ||
136 | CHECK_ABSOLUTE_ERROR(f2(x), cheb2(x), acceptable_error_2); | |
b32b8144 FG |
137 | x += static_cast<Real>(1)/static_cast<Real>(1 << 7); |
138 | } | |
139 | } | |
140 | ||
20effc67 | 141 | |
b32b8144 FG |
142 | //Validate that the Chebyshev polynomials are well approximated by the Chebyshev transform. |
143 | template<class Real> | |
144 | void test_chebyshev_chebyshev_transform() | |
145 | { | |
146 | Real tol = 500*std::numeric_limits<Real>::epsilon(); | |
147 | // T_0 = 1: | |
148 | auto t0 = [](Real) { return 1; }; | |
149 | chebyshev_transform<Real> cheb0(t0, -1, 1); | |
20effc67 | 150 | CHECK_ABSOLUTE_ERROR(2, cheb0.coefficients()[0], tol); |
b32b8144 FG |
151 | |
152 | Real x = -1; | |
153 | while (x < 1) | |
154 | { | |
20effc67 TL |
155 | CHECK_ABSOLUTE_ERROR(1, cheb0(x), tol); |
156 | CHECK_ABSOLUTE_ERROR(Real(0), cheb0.prime(x), tol); | |
b32b8144 FG |
157 | x += static_cast<Real>(1)/static_cast<Real>(1 << 7); |
158 | } | |
159 | ||
160 | // T_1 = x: | |
161 | auto t1 = [](Real x) { return x; }; | |
162 | chebyshev_transform<Real> cheb1(t1, -1, 1); | |
20effc67 | 163 | CHECK_ABSOLUTE_ERROR(Real(1), cheb1.coefficients()[1], tol); |
b32b8144 FG |
164 | |
165 | x = -1; | |
166 | while (x < 1) | |
167 | { | |
20effc67 TL |
168 | CHECK_ABSOLUTE_ERROR(x, cheb1(x), tol); |
169 | CHECK_ABSOLUTE_ERROR(Real(1), cheb1.prime(x), tol); | |
b32b8144 FG |
170 | x += static_cast<Real>(1)/static_cast<Real>(1 << 7); |
171 | } | |
172 | ||
173 | ||
174 | auto t2 = [](Real x) { return 2*x*x-1; }; | |
175 | chebyshev_transform<Real> cheb2(t2, -1, 1); | |
20effc67 | 176 | CHECK_ABSOLUTE_ERROR(Real(1), cheb2.coefficients()[2], tol); |
b32b8144 FG |
177 | |
178 | x = -1; | |
179 | while (x < 1) | |
180 | { | |
20effc67 TL |
181 | CHECK_ABSOLUTE_ERROR(t2(x), cheb2(x), tol); |
182 | CHECK_ABSOLUTE_ERROR(4*x, cheb2.prime(x), tol); | |
b32b8144 FG |
183 | x += static_cast<Real>(1)/static_cast<Real>(1 << 7); |
184 | } | |
185 | } | |
186 | ||
20effc67 | 187 | int main() |
b32b8144 FG |
188 | { |
189 | #ifdef TEST1 | |
190 | test_chebyshev_chebyshev_transform<float>(); | |
191 | test_sin_chebyshev_transform<float>(); | |
192 | test_atap_examples<float>(); | |
193 | test_sinc_chebyshev_transform<float>(); | |
194 | #endif | |
195 | #ifdef TEST2 | |
196 | test_chebyshev_chebyshev_transform<double>(); | |
197 | test_sin_chebyshev_transform<double>(); | |
198 | test_atap_examples<double>(); | |
199 | test_sinc_chebyshev_transform<double>(); | |
200 | #endif | |
201 | #ifdef TEST3 | |
202 | test_chebyshev_chebyshev_transform<long double>(); | |
203 | test_sin_chebyshev_transform<long double>(); | |
204 | test_atap_examples<long double>(); | |
205 | test_sinc_chebyshev_transform<long double>(); | |
206 | #endif | |
207 | #ifdef TEST4 | |
208 | #ifdef BOOST_HAS_FLOAT128 | |
209 | test_chebyshev_chebyshev_transform<__float128>(); | |
210 | test_sin_chebyshev_transform<__float128>(); | |
211 | test_atap_examples<__float128>(); | |
212 | test_sinc_chebyshev_transform<__float128>(); | |
213 | #endif | |
214 | #endif | |
20effc67 TL |
215 | |
216 | return boost::math::test::report_errors(); | |
b32b8144 FG |
217 | } |
218 | ||
219 |