]>
Commit | Line | Data |
---|---|---|
11fdf7f2 TL |
1 | // (C) Copyright Nick Thompson, 2018 |
2 | // Use, modification and distribution are subject to the | |
3 | // Boost Software License, Version 1.0. (See accompanying file | |
4 | // LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt) | |
5 | ||
6 | #define BOOST_TEST_MODULE numerical_differentiation_test | |
7 | ||
8 | #include <cmath> | |
9 | #include <limits> | |
10 | #include <iostream> | |
11 | #include <boost/type_index.hpp> | |
12 | #include <boost/test/included/unit_test.hpp> | |
92f5a8d4 | 13 | #include <boost/test/tools/floating_point_comparison.hpp> |
11fdf7f2 TL |
14 | #include <boost/math/special_functions/bessel.hpp> |
15 | #include <boost/math/special_functions/bessel_prime.hpp> | |
16 | #include <boost/math/special_functions/next.hpp> | |
92f5a8d4 | 17 | #include <boost/math/differentiation/finite_difference.hpp> |
11fdf7f2 TL |
18 | |
19 | using std::abs; | |
20 | using std::pow; | |
92f5a8d4 TL |
21 | using boost::math::differentiation::finite_difference_derivative; |
22 | using boost::math::differentiation::complex_step_derivative; | |
11fdf7f2 TL |
23 | using boost::math::cyl_bessel_j; |
24 | using boost::math::cyl_bessel_j_prime; | |
25 | using boost::math::constants::half; | |
11fdf7f2 TL |
26 | |
27 | template<class Real, size_t order> | |
28 | void test_order(size_t points_to_test) | |
29 | { | |
30 | std::cout << "Testing order " << order << " derivative error estimate on type " << boost::typeindex::type_id<Real>().pretty_name() << "\n"; | |
31 | std::cout << std::setprecision(std::numeric_limits<Real>::digits10); | |
32 | //std::cout << std::fixed << std::scientific; | |
92f5a8d4 | 33 | auto f = [](Real t) { return boost::math::cyl_bessel_j<Real>(1, t); }; |
11fdf7f2 TL |
34 | Real min = -100000.0; |
35 | Real max = -min; | |
36 | Real x = min; | |
37 | Real max_error = 0; | |
38 | Real max_relative_error_in_error = 0; | |
39 | size_t j = 0; | |
40 | size_t failures = 0; | |
41 | while (j < points_to_test) | |
42 | { | |
43 | x = min + (Real) 2*j*max/ (Real) points_to_test; | |
44 | Real error_estimate; | |
45 | Real computed = finite_difference_derivative<decltype(f), Real, order>(f, x, &error_estimate); | |
46 | Real expected = (Real) cyl_bessel_j_prime<Real>(1, x); | |
47 | Real error = abs(computed - expected); | |
48 | // The error estimate is provided under the assumption that the function is evaluated to 1 ULP. | |
49 | // Presumably no one will be too offended by this estimate being off by a factor of 2 or so. | |
50 | if (error > 2*error_estimate) | |
51 | { | |
52 | ++failures; | |
53 | Real relative_error_in_error = abs(error - error_estimate)/ error; | |
54 | if (relative_error_in_error > max_relative_error_in_error) | |
55 | { | |
56 | max_relative_error_in_error = relative_error_in_error; | |
57 | } | |
58 | if (relative_error_in_error > 2) | |
59 | { | |
60 | throw std::logic_error("Relative error in error is too high!"); | |
61 | } | |
62 | } | |
63 | if (error > max_error) | |
64 | { | |
65 | max_error = error; | |
66 | } | |
67 | ++j; | |
68 | } | |
69 | //std::cout << "Maximum error :" << max_error << "\n"; | |
70 | //std::cout << "Error estimate failed " << failures << " times out of " << points_to_test << "\n"; | |
71 | //std::cout << "Failure rate: " << (double) failures / (double) points_to_test << "\n"; | |
72 | //std::cout << "Maximum error in estimated error = " << max_relative_error_in_error << "\n"; | |
73 | //Real convergence_rate = (Real) order/ (Real) (order + 1); | |
74 | //std::cout << "eps^(order/order+1) = " << pow(std::numeric_limits<Real>::epsilon(), convergence_rate) << "\n\n\n"; | |
75 | ||
76 | bool max_error_good = max_error < 2*sqrt(std::numeric_limits<Real>::epsilon()); | |
77 | BOOST_TEST(max_error_good); | |
78 | ||
79 | bool error_estimate_good = max_relative_error_in_error < (Real) 2; | |
80 | BOOST_TEST(error_estimate_good); | |
81 | ||
82 | double failure_rate = (double) failures / (double) points_to_test; | |
83 | BOOST_CHECK_SMALL(failure_rate, 0.05); | |
84 | } | |
85 | ||
86 | template<class Real> | |
87 | void test_bessel() | |
88 | { | |
89 | std::cout << "Testing numerical derivatives of Bessel's function on type " << boost::typeindex::type_id<Real>().pretty_name() << "\n"; | |
90 | std::cout << std::setprecision(std::numeric_limits<Real>::digits10); | |
91 | ||
92 | Real eps = std::numeric_limits<Real>::epsilon(); | |
93 | Real x = static_cast<Real>(25.1); | |
92f5a8d4 | 94 | auto f = [](Real t) { return boost::math::cyl_bessel_j(12, t); }; |
11fdf7f2 TL |
95 | |
96 | Real computed = finite_difference_derivative<decltype(f), Real, 1>(f, x); | |
97 | Real expected = cyl_bessel_j_prime(12, x); | |
98 | Real error_estimate = 4*abs(f(x))*sqrt(eps); | |
99 | //std::cout << std::setprecision(std::numeric_limits<Real>::digits10); | |
100 | //std::cout << "cyl_bessel_j_prime: " << expected << std::endl; | |
101 | //std::cout << "First order fd : " << computed << std::endl; | |
102 | //std::cout << "Error : " << abs(computed - expected) << std::endl; | |
103 | //std::cout << "a prior error est : " << error_estimate << std::endl; | |
104 | ||
105 | BOOST_CHECK_CLOSE_FRACTION(expected, computed, 10*error_estimate); | |
106 | ||
107 | computed = finite_difference_derivative<decltype(f), Real, 2>(f, x); | |
108 | expected = cyl_bessel_j_prime(12, x); | |
109 | error_estimate = abs(f(x))*pow(eps, boost::math::constants::two_thirds<Real>()); | |
110 | //std::cout << std::setprecision(std::numeric_limits<Real>::digits10); | |
111 | //std::cout << "cyl_bessel_j_prime: " << expected << std::endl; | |
112 | //std::cout << "Second order fd : " << computed << std::endl; | |
113 | //std::cout << "Error : " << abs(computed - expected) << std::endl; | |
114 | //std::cout << "a prior error est : " << error_estimate << std::endl; | |
115 | ||
116 | BOOST_CHECK_CLOSE_FRACTION(expected, computed, 50*error_estimate); | |
117 | ||
118 | computed = finite_difference_derivative<decltype(f), Real, 4>(f, x); | |
119 | expected = cyl_bessel_j_prime(12, x); | |
120 | error_estimate = abs(f(x))*pow(eps, (Real) 4 / (Real) 5); | |
121 | //std::cout << std::setprecision(std::numeric_limits<Real>::digits10); | |
122 | //std::cout << "cyl_bessel_j_prime: " << expected << std::endl; | |
123 | //std::cout << "Fourth order fd : " << computed << std::endl; | |
124 | //std::cout << "Error : " << abs(computed - expected) << std::endl; | |
125 | //std::cout << "a prior error est : " << error_estimate << std::endl; | |
126 | ||
127 | BOOST_CHECK_CLOSE_FRACTION(expected, computed, 25*error_estimate); | |
128 | ||
129 | ||
130 | computed = finite_difference_derivative<decltype(f), Real, 6>(f, x); | |
131 | expected = cyl_bessel_j_prime(12, x); | |
132 | error_estimate = abs(f(x))*pow(eps, (Real) 6/ (Real) 7); | |
133 | //std::cout << std::setprecision(std::numeric_limits<Real>::digits10); | |
134 | //std::cout << "cyl_bessel_j_prime: " << expected << std::endl; | |
135 | //std::cout << "Sixth order fd : " << computed << std::endl; | |
136 | //std::cout << "Error : " << abs(computed - expected) << std::endl; | |
137 | //std::cout << "a prior error est : " << error_estimate << std::endl; | |
138 | ||
139 | BOOST_CHECK_CLOSE_FRACTION(expected, computed, 100*error_estimate); | |
140 | ||
141 | computed = finite_difference_derivative<decltype(f), Real, 8>(f, x); | |
142 | expected = cyl_bessel_j_prime(12, x); | |
143 | error_estimate = abs(f(x))*pow(eps, (Real) 8/ (Real) 9); | |
144 | //std::cout << std::setprecision(std::numeric_limits<Real>::digits10); | |
145 | //std::cout << "cyl_bessel_j_prime: " << expected << std::endl; | |
146 | //std::cout << "Eighth order fd : " << computed << std::endl; | |
147 | //std::cout << "Error : " << abs(computed - expected) << std::endl; | |
148 | //std::cout << "a prior error est : " << error_estimate << std::endl; | |
149 | ||
150 | BOOST_CHECK_CLOSE_FRACTION(expected, computed, 25*error_estimate); | |
151 | } | |
152 | ||
153 | // Example of a function which is subject to catastrophic cancellation using finite-differences, but is almost perfectly stable using complex step: | |
154 | template<class RealOrComplex> | |
155 | RealOrComplex moler_example(RealOrComplex x) | |
156 | { | |
157 | using std::sin; | |
158 | using std::cos; | |
159 | using std::exp; | |
160 | ||
161 | RealOrComplex cosx = cos(x); | |
162 | RealOrComplex sinx = sin(x); | |
163 | return exp(x)/(cosx*cosx*cosx + sinx*sinx*sinx); | |
164 | } | |
165 | ||
166 | template<class RealOrComplex> | |
167 | RealOrComplex moler_example_derivative(RealOrComplex x) | |
168 | { | |
169 | using std::sin; | |
170 | using std::cos; | |
171 | using std::exp; | |
172 | ||
173 | RealOrComplex expx = exp(x); | |
174 | RealOrComplex cosx = cos(x); | |
175 | RealOrComplex sinx = sin(x); | |
176 | RealOrComplex coscubed_sincubed = cosx*cosx*cosx + sinx*sinx*sinx; | |
177 | return (expx/coscubed_sincubed)*(1 - 3*(sinx*sinx*cosx - sinx*cosx*cosx)/ (coscubed_sincubed)); | |
178 | } | |
179 | ||
180 | ||
181 | template<class Real> | |
182 | void test_complex_step() | |
183 | { | |
184 | using std::abs; | |
185 | using std::complex; | |
186 | using std::isfinite; | |
187 | using std::isnormal; | |
188 | std::cout << "Testing numerical derivatives of Bessel's function on type " << boost::typeindex::type_id<Real>().pretty_name() << "\n"; | |
189 | std::cout << std::setprecision(std::numeric_limits<Real>::digits10); | |
190 | Real x = -100; | |
191 | while ( x < 100 ) | |
192 | { | |
193 | if (!isfinite(moler_example(x))) | |
194 | { | |
195 | x += 1; | |
196 | continue; | |
197 | } | |
198 | Real expected = moler_example_derivative<Real>(x); | |
199 | Real computed = complex_step_derivative(moler_example<complex<Real>>, x); | |
200 | if (!isfinite(expected)) | |
201 | { | |
202 | x += 1; | |
203 | continue; | |
204 | } | |
205 | if (abs(expected) <= std::numeric_limits<Real>::epsilon()) | |
206 | { | |
207 | bool issmall = computed < std::numeric_limits<Real>::epsilon(); | |
208 | BOOST_TEST(issmall); | |
209 | } | |
210 | else | |
211 | { | |
212 | BOOST_CHECK_CLOSE_FRACTION(expected, computed, 200*std::numeric_limits<Real>::epsilon()); | |
213 | } | |
214 | x += 1; | |
215 | } | |
216 | } | |
217 | ||
218 | ||
219 | BOOST_AUTO_TEST_CASE(numerical_differentiation_test) | |
220 | { | |
221 | test_complex_step<float>(); | |
222 | test_complex_step<double>(); | |
92f5a8d4 | 223 | |
11fdf7f2 TL |
224 | test_bessel<float>(); |
225 | test_bessel<double>(); | |
92f5a8d4 | 226 | |
11fdf7f2 TL |
227 | |
228 | size_t points_to_test = 1000; | |
229 | test_order<float, 1>(points_to_test); | |
230 | test_order<double, 1>(points_to_test); | |
92f5a8d4 | 231 | |
11fdf7f2 TL |
232 | |
233 | test_order<float, 2>(points_to_test); | |
234 | test_order<double, 2>(points_to_test); | |
11fdf7f2 TL |
235 | |
236 | test_order<float, 4>(points_to_test); | |
237 | test_order<double, 4>(points_to_test); | |
11fdf7f2 TL |
238 | |
239 | test_order<float, 6>(points_to_test); | |
240 | test_order<double, 6>(points_to_test); | |
11fdf7f2 TL |
241 | |
242 | test_order<float, 8>(points_to_test); | |
243 | test_order<double, 8>(points_to_test); | |
92f5a8d4 | 244 | |
11fdf7f2 | 245 | } |