boost_mirror/math
1
0
Fork 0
mirror of https://github.com/boostorg/math.git synced 2025-05-11 21:33:52 +00:00
Code Issues Packages Projects Releases Wiki Activity

Add test and cast to finite differences

Browse Source
This commit is contained in:
Matt Borland 2023-05-30 16:46:23 +02:00
parent d1bd7b6f80
commit 823fcd4cf2
No known key found for this signature in database
GPG Key ID: 1FC9FE1989F22A0B
2 changed files with 37 additions and 20 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

4
include/boost/math/differentiation/finite_difference.hpp
View File

@ -153,7 +153,7 @@ namespace detail {
const Real eps = (numeric_limits<Real>::epsilon)();
// Error bound ~eps^4/5
Real h = pow(11.25*eps, static_cast<Real>(1) / static_cast<Real>(5));
Real h = pow(Real(11.25)*eps, static_cast<Real>(1) / static_cast<Real>(5));
h = detail::make_xph_representable(x, h);
Real ymth = f(x - 2 * h);
Real yth = f(x + 2 * h);
@ -222,7 +222,7 @@ namespace detail {
// Mathematica code to get the error:
// Series[(f[x+h]-f[x-h])*(4/5) + (1/5)*(f[x-2*h] - f[x+2*h]) + (4/105)*(f[x+3*h] - f[x-3*h]) + (1/280)*(f[x-4*h] - f[x+4*h]), {h, 0, 9}]
// If we used Kahan summation, we could get the max error down to h^8|f^(9)(x)|/630 + |f(x)|eps/h.
Real h = pow(551.25*eps, static_cast<Real>(1) / static_cast<Real>(9));
Real h = pow(Real(551.25)*eps, static_cast<Real>(1) / static_cast<Real>(9));
h = detail::make_xph_representable(x, h);
Real yh = f(x + h);

53
test/test_numerical_differentiation.cpp
View File

@ -16,6 +16,10 @@
#include <boost/math/special_functions/next.hpp>
#include <boost/math/differentiation/finite_difference.hpp>
#if __has_include(<stdfloat>)
# include <stdfloat>
#endif
using std::abs;
using std::pow;
using boost::math::differentiation::finite_difference_derivative;
@ -31,7 +35,7 @@ void test_order(size_t points_to_test)
std::cout << std::setprecision(std::numeric_limits<Real>::digits10);
//std::cout << std::fixed << std::scientific;
auto f = [](Real t) { return boost::math::cyl_bessel_j<Real>(1, t); };
Real min = -100000.0;
Real min = Real(-100000.0);
Real max = -min;
Real x = min;
Real max_error = 0;
@ -95,7 +99,7 @@ void test_bessel()
Real computed = finite_difference_derivative<decltype(f), Real, 1>(f, x);
Real expected = cyl_bessel_j_prime(12, x);
Real error_estimate = 4*abs(f(x))*sqrt(eps);
Real error_estimate = Real(4*abs(f(x))*sqrt(eps));
//std::cout << std::setprecision(std::numeric_limits<Real>::digits10);
//std::cout << "cyl_bessel_j_prime: " << expected << std::endl;
//std::cout << "First order fd : " << computed << std::endl;
@ -218,28 +222,41 @@ void test_complex_step()
BOOST_AUTO_TEST_CASE(numerical_differentiation_test)
{
constexpr size_t points_to_test = 1000;
#ifdef __STDCPP_FLOAT32_T__
test_complex_step<std::float32_t>();
test_bessel<std::float32_t>();
test_order<std::float32_t, 1>(points_to_test);
test_order<std::float32_t, 2>(points_to_test);
test_order<std::float32_t, 4>(points_to_test);
test_order<std::float32_t, 6>(points_to_test);
test_order<std::float32_t, 8>(points_to_test);
#else
test_complex_step<float>();
test_complex_step<double>();
test_bessel<float>();
test_bessel<double>();
size_t points_to_test = 1000;
test_order<float, 1>(points_to_test);
test_order<double, 1>(points_to_test);
test_order<float, 2>(points_to_test);
test_order<double, 2>(points_to_test);
test_order<float, 4>(points_to_test);
test_order<double, 4>(points_to_test);
test_order<float, 6>(points_to_test);
test_order<double, 6>(points_to_test);
test_order<float, 8>(points_to_test);
test_order<double, 8>(points_to_test);
#endif
#ifdef __STDCPP_FLOAT64_T__
test_complex_step<std::float64_t>();
test_bessel<std::float64_t>();
test_order<std::float64_t, 1>(points_to_test);
test_order<std::float64_t, 2>(points_to_test);
test_order<std::float64_t, 4>(points_to_test);
test_order<std::float64_t, 6>(points_to_test);
test_order<std::float64_t, 8>(points_to_test);
#else
test_complex_step<double>();
test_bessel<double>();
test_order<double, 1>(points_to_test);
test_order<double, 2>(points_to_test);
test_order<double, 4>(points_to_test);
test_order<double, 6>(points_to_test);
test_order<double, 8>(points_to_test);
#endif
}