libs/math/test/test_numerical_differentiation.cpp

   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>
  13 #include <boost/test/tools/floating_point_comparison.hpp>
  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>
  17 #include <boost/math/differentiation/finite_difference.hpp>
  18
  19 using std::abs;
  20 using std::pow;
  21 using boost::math::differentiation::finite_difference_derivative;
  22 using boost::math::differentiation::complex_step_derivative;
  23 using boost::math::cyl_bessel_j;
  24 using boost::math::cyl_bessel_j_prime;
  25 using boost::math::constants::half;
  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;
  33     auto f = [](Real t) { return boost::math::cyl_bessel_j<Real>(1, t); };
  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);
  94     auto f = [](Real t) { return boost::math::cyl_bessel_j(12, t); };
  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>();
 223
 224     test_bessel<float>();
 225     test_bessel<double>();
 226
 227
 228     size_t points_to_test = 1000;
 229     test_order<float, 1>(points_to_test);
 230     test_order<double, 1>(points_to_test);
 231
 232
 233     test_order<float, 2>(points_to_test);
 234     test_order<double, 2>(points_to_test);
 235
 236     test_order<float, 4>(points_to_test);
 237     test_order<double, 4>(points_to_test);
 238
 239     test_order<float, 6>(points_to_test);
 240     test_order<double, 6>(points_to_test);
 241
 242     test_order<float, 8>(points_to_test);
 243     test_order<double, 8>(points_to_test);
 244
 245 }