267 lines · plain
1// (C) Copyright Nick Thompson 2018.2// Use, modification and distribution are subject to the3// Boost Software License, Version 1.0. (See accompanying file4// LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt)5 6#ifndef BOOST_MATH_DIFFERENTIATION_FINITE_DIFFERENCE_HPP7#define BOOST_MATH_DIFFERENTIATION_FINITE_DIFFERENCE_HPP8 9/*10 * Performs numerical differentiation by finite-differences.11 *12 * All numerical differentiation using finite-differences are ill-conditioned, and these routines are no exception.13 * A simple argument demonstrates that the error is unbounded as h->0.14 * Take the one sides finite difference formula f'(x) = (f(x+h)-f(x))/h.15 * The evaluation of f induces an error as well as the error from the finite-difference approximation, giving16 * |f'(x) - (f(x+h) -f(x))/h| < h|f''(x)|/2 + (|f(x)|+|f(x+h)|)eps/h =: g(h), where eps is the unit roundoff for the type.17 * It is reasonable to choose h in a way that minimizes the maximum error bound g(h).18 * The value of h that minimizes g is h = sqrt(2eps(|f(x)| + |f(x+h)|)/|f''(x)|), and for this value of h the error bound is19 * sqrt(2eps(|f(x+h) +f(x)||f''(x)|)).20 * In fact it is not necessary to compute the ratio (|f(x+h)| + |f(x)|)/|f''(x)|; the error bound of ~\sqrt{\epsilon} still holds if we set it to one.21 *22 *23 * For more details on this method of analysis, see24 *25 * http://www.uio.no/studier/emner/matnat/math/MAT-INF1100/h08/kompendiet/diffint.pdf26 * http://web.archive.org/web/20150420195907/http://www.uio.no/studier/emner/matnat/math/MAT-INF1100/h08/kompendiet/diffint.pdf27 *28 *29 * It can be shown on general grounds that when choosing the optimal h, the maximum error in f'(x) is ~(|f(x)|eps)^k/k+1|f^(k-1)(x)|^1/k+1.30 * From this we can see that full precision can be recovered in the limit k->infinity.31 *32 * References:33 *34 * 1) Fornberg, Bengt. "Generation of finite difference formulas on arbitrarily spaced grids." Mathematics of computation 51.184 (1988): 699-706.35 *36 *37 * The second algorithm, the complex step derivative, is not ill-conditioned.38 * However, it requires that your function can be evaluated at complex arguments.39 * The idea is that f(x+ih) = f(x) +ihf'(x) - h^2f''(x) + ... so f'(x) \approx Im[f(x+ih)]/h.40 * No subtractive cancellation occurs. The error is ~ eps|f'(x)| + eps^2|f'''(x)|/6; hard to beat that.41 *42 * References:43 *44 * 1) Squire, William, and George Trapp. "Using complex variables to estimate derivatives of real functions." Siam Review 40.1 (1998): 110-112.45 */46 47#include <complex>48#include <boost/math/special_functions/next.hpp>49 50namespace boost{ namespace math{ namespace differentiation {51 52namespace detail {53 template<class Real>54 Real make_xph_representable(Real x, Real h)55 {56 using std::numeric_limits;57 // Redefine h so that x + h is representable. Not using this trick leads to large error.58 // The compiler flag -ffast-math evaporates these operations . . .59 Real temp = x + h;60 h = temp - x;61 // Handle the case x + h == x:62 if (h == 0)63 {64 h = boost::math::nextafter(x, (numeric_limits<Real>::max)()) - x;65 }66 return h;67 }68}69 70template<class F, class Real>71Real complex_step_derivative(const F f, Real x)72{73 // Is it really this easy? Yes.74 // Note that some authors recommend taking the stepsize h to be smaller than epsilon(), some recommending use of the min().75 // This idea was tested over a few billion test cases and found the make the error *much* worse.76 // Even 2eps and eps/2 made the error worse, which was surprising.77 using std::complex;78 using std::numeric_limits;79 constexpr const Real step = (numeric_limits<Real>::epsilon)();80 constexpr const Real inv_step = 1/(numeric_limits<Real>::epsilon)();81 return f(complex<Real>(x, step)).imag()*inv_step;82}83 84namespace detail {85 86 template <unsigned>87 struct fd_tag {};88 89 template<class F, class Real>90 Real finite_difference_derivative(const F f, Real x, Real* error, const fd_tag<1>&)91 {92 using std::sqrt;93 using std::pow;94 using std::abs;95 using std::numeric_limits;96 97 const Real eps = (numeric_limits<Real>::epsilon)();98 // Error bound ~eps^1/299 // Note that this estimate of h differs from the best estimate by a factor of sqrt((|f(x)| + |f(x+h)|)/|f''(x)|).100 // Since this factor is invariant under the scaling f -> kf, then we are somewhat justified in approximating it by 1.101 // This approximation will get better as we move to higher orders of accuracy.102 Real h = 2 * sqrt(eps);103 h = detail::make_xph_representable(x, h);104 105 Real yh = f(x + h);106 Real y0 = f(x);107 Real diff = yh - y0;108 if (error)109 {110 Real ym = f(x - h);111 Real ypph = abs(yh - 2 * y0 + ym) / h;112 // h*|f''(x)|*0.5 + (|f(x+h)+|f(x)|)*eps/h113 *error = ypph / 2 + (abs(yh) + abs(y0))*eps / h;114 }115 return diff / h;116 }117 118 template<class F, class Real>119 Real finite_difference_derivative(const F f, Real x, Real* error, const fd_tag<2>&)120 {121 using std::sqrt;122 using std::pow;123 using std::abs;124 using std::numeric_limits;125 126 const Real eps = (numeric_limits<Real>::epsilon)();127 // Error bound ~eps^2/3128 // See the previous discussion to understand determination of h and the error bound.129 // Series[(f[x+h] - f[x-h])/(2*h), {h, 0, 4}]130 Real h = pow(3 * eps, static_cast<Real>(1) / static_cast<Real>(3));131 h = detail::make_xph_representable(x, h);132 133 Real yh = f(x + h);134 Real ymh = f(x - h);135 Real diff = yh - ymh;136 if (error)137 {138 Real yth = f(x + 2 * h);139 Real ymth = f(x - 2 * h);140 *error = eps * (abs(yh) + abs(ymh)) / (2 * h) + abs((yth - ymth) / 2 - diff) / (6 * h);141 }142 143 return diff / (2 * h);144 }145 146 template<class F, class Real>147 Real finite_difference_derivative(const F f, Real x, Real* error, const fd_tag<4>&)148 {149 using std::sqrt;150 using std::pow;151 using std::abs;152 using std::numeric_limits;153 154 const Real eps = (numeric_limits<Real>::epsilon)();155 // Error bound ~eps^4/5156 Real h = pow(Real(11.25)*eps, static_cast<Real>(1) / static_cast<Real>(5));157 h = detail::make_xph_representable(x, h);158 Real ymth = f(x - 2 * h);159 Real yth = f(x + 2 * h);160 Real yh = f(x + h);161 Real ymh = f(x - h);162 Real y2 = ymth - yth;163 Real y1 = yh - ymh;164 if (error)165 {166 // Mathematica code to extract the remainder:167 // Series[(f[x-2*h]+ 8*f[x+h] - 8*f[x-h] - f[x+2*h])/(12*h), {h, 0, 7}]168 Real y_three_h = f(x + 3 * h);169 Real y_m_three_h = f(x - 3 * h);170 // Error from fifth derivative:171 *error = abs((y_three_h - y_m_three_h) / 2 + 2 * (ymth - yth) + 5 * (yh - ymh) / 2) / (30 * h);172 // Error from function evaluation:173 *error += eps * (abs(yth) + abs(ymth) + 8 * (abs(ymh) + abs(yh))) / (12 * h);174 }175 return (y2 + 8 * y1) / (12 * h);176 }177 178 template<class F, class Real>179 Real finite_difference_derivative(const F f, Real x, Real* error, const fd_tag<6>&)180 {181 using std::sqrt;182 using std::pow;183 using std::abs;184 using std::numeric_limits;185 186 const Real eps = (numeric_limits<Real>::epsilon)();187 // Error bound ~eps^6/7188 // Error: h^6f^(7)(x)/140 + 5|f(x)|eps/h189 Real h = pow(eps / 168, static_cast<Real>(1) / static_cast<Real>(7));190 h = detail::make_xph_representable(x, h);191 192 Real yh = f(x + h);193 Real ymh = f(x - h);194 Real y1 = yh - ymh;195 Real y2 = f(x - 2 * h) - f(x + 2 * h);196 Real y3 = f(x + 3 * h) - f(x - 3 * h);197 198 if (error)199 {200 // Mathematica code to generate fd scheme for 7th derivative:201 // Sum[(-1)^i*Binomial[7, i]*(f[x+(3-i)*h] + f[x+(4-i)*h])/2, {i, 0, 7}]202 // Mathematica to demonstrate that this is a finite difference formula for 7th derivative:203 // Series[(f[x+4*h]-f[x-4*h] + 6*(f[x-3*h] - f[x+3*h]) + 14*(f[x-h] - f[x+h] + f[x+2*h] - f[x-2*h]))/2, {h, 0, 15}]204 Real y7 = (f(x + 4 * h) - f(x - 4 * h) - 6 * y3 - 14 * y1 - 14 * y2) / 2;205 *error = abs(y7) / (140 * h) + 5 * (abs(yh) + abs(ymh))*eps / h;206 }207 return (y3 + 9 * y2 + 45 * y1) / (60 * h);208 }209 210 template<class F, class Real>211 Real finite_difference_derivative(const F f, Real x, Real* error, const fd_tag<8>&)212 {213 using std::sqrt;214 using std::pow;215 using std::abs;216 using std::numeric_limits;217 218 const Real eps = (numeric_limits<Real>::epsilon)();219 // Error bound ~eps^8/9.220 // In double precision, we only expect to lose two digits of precision while using this formula, at the cost of 8 function evaluations.221 // Error: h^8|f^(9)(x)|/630 + 7|f(x)|eps/h assuming 7 unstabilized additions.222 // Mathematica code to get the error:223 // 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}]224 // If we used Kahan summation, we could get the max error down to h^8|f^(9)(x)|/630 + |f(x)|eps/h.225 Real h = pow(Real(551.25)*eps, static_cast<Real>(1) / static_cast<Real>(9));226 h = detail::make_xph_representable(x, h);227 228 Real yh = f(x + h);229 Real ymh = f(x - h);230 Real y1 = yh - ymh;231 Real y2 = f(x - 2 * h) - f(x + 2 * h);232 Real y3 = f(x + 3 * h) - f(x - 3 * h);233 Real y4 = f(x - 4 * h) - f(x + 4 * h);234 235 Real tmp1 = 3 * y4 / 8 + 4 * y3;236 Real tmp2 = 21 * y2 + 84 * y1;237 238 if (error)239 {240 // Mathematica code to generate fd scheme for 7th derivative:241 // Sum[(-1)^i*Binomial[9, i]*(f[x+(4-i)*h] + f[x+(5-i)*h])/2, {i, 0, 9}]242 // Mathematica to demonstrate that this is a finite difference formula for 7th derivative:243 // Series[(f[x+5*h]-f[x- 5*h])/2 + 4*(f[x-4*h] - f[x+4*h]) + 27*(f[x+3*h] - f[x-3*h])/2 + 24*(f[x-2*h] - f[x+2*h]) + 21*(f[x+h] - f[x-h]), {h, 0, 15}]244 Real f9 = (f(x + 5 * h) - f(x - 5 * h)) / 2 + 4 * y4 + 27 * y3 / 2 + 24 * y2 + 21 * y1;245 *error = abs(f9) / (630 * h) + 7 * (abs(yh) + abs(ymh))*eps / h;246 }247 return (tmp1 + tmp2) / (105 * h);248 }249 250 template<class F, class Real, class tag>251 Real finite_difference_derivative(const F, Real, Real*, const tag&)252 {253 // Always fails, but condition is template-arg-dependent so only evaluated if we get instantiated.254 static_assert(sizeof(Real) == 0, "Finite difference not implemented for this order: try 1, 2, 4, 6 or 8");255 }256 257}258 259template<class F, class Real, size_t order=6>260inline Real finite_difference_derivative(const F f, Real x, Real* error = nullptr)261{262 return detail::finite_difference_derivative(f, x, error, detail::fd_tag<order>());263}264 265}}} // namespaces266#endif267