math/differentiation/finite_difference.hpp

0001 //  (C) Copyright Nick Thompson 2018.
0002 //  Use, modification and distribution are subject to the
0003 //  Boost Software License, Version 1.0. (See accompanying file
0004 //  LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt)
0005
0006 #ifndef BOOST_MATH_DIFFERENTIATION_FINITE_DIFFERENCE_HPP
0007 #define BOOST_MATH_DIFFERENTIATION_FINITE_DIFFERENCE_HPP
0008
0009 /*
0010  * Performs numerical differentiation by finite-differences.
0011  *
0012  * All numerical differentiation using finite-differences are ill-conditioned, and these routines are no exception.
0013  * A simple argument demonstrates that the error is unbounded as h->0.
0014  * Take the one sides finite difference formula f'(x) = (f(x+h)-f(x))/h.
0015  * The evaluation of f induces an error as well as the error from the finite-difference approximation, giving
0016  * |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.
0017  * It is reasonable to choose h in a way that minimizes the maximum error bound g(h).
0018  * 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 is
0019  * sqrt(2eps(|f(x+h) +f(x)||f''(x)|)).
0020  * 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.
0021  *
0022  *
0023  * For more details on this method of analysis, see
0024  *
0025  * http://www.uio.no/studier/emner/matnat/math/MAT-INF1100/h08/kompendiet/diffint.pdf
0026  * http://web.archive.org/web/20150420195907/http://www.uio.no/studier/emner/matnat/math/MAT-INF1100/h08/kompendiet/diffint.pdf
0027  *
0028  *
0029  * 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.
0030  * From this we can see that full precision can be recovered in the limit k->infinity.
0031  *
0032  * References:
0033  *
0034  * 1) Fornberg, Bengt. "Generation of finite difference formulas on arbitrarily spaced grids." Mathematics of computation 51.184 (1988): 699-706.
0035  *
0036  *
0037  * The second algorithm, the complex step derivative, is not ill-conditioned.
0038  * However, it requires that your function can be evaluated at complex arguments.
0039  * The idea is that f(x+ih) = f(x) +ihf'(x) - h^2f''(x) + ... so f'(x) \approx Im[f(x+ih)]/h.
0040  * No subtractive cancellation occurs. The error is ~ eps|f'(x)| + eps^2|f'''(x)|/6; hard to beat that.
0041  *
0042  * References:
0043  *
0044  * 1) Squire, William, and George Trapp. "Using complex variables to estimate derivatives of real functions." Siam Review 40.1 (1998): 110-112.
0045  */
0046
0047 #include <complex>
0048 #include <boost/math/special_functions/next.hpp>
0049
0050 namespace boost{ namespace math{ namespace differentiation {
0051
0052 namespace detail {
0053     template<class Real>
0054     Real make_xph_representable(Real x, Real h)
0055     {
0056         using std::numeric_limits;
0057         // Redefine h so that x + h is representable. Not using this trick leads to large error.
0058         // The compiler flag -ffast-math evaporates these operations . . .
0059         Real temp = x + h;
0060         h = temp - x;
0061         // Handle the case x + h == x:
0062         if (h == 0)
0063         {
0064             h = boost::math::nextafter(x, (numeric_limits<Real>::max)()) - x;
0065         }
0066         return h;
0067     }
0068 }
0069
0070 template<class F, class Real>
0071 Real complex_step_derivative(const F f, Real x)
0072 {
0073     // Is it really this easy? Yes.
0074     // Note that some authors recommend taking the stepsize h to be smaller than epsilon(), some recommending use of the min().
0075     // This idea was tested over a few billion test cases and found the make the error *much* worse.
0076     // Even 2eps and eps/2 made the error worse, which was surprising.
0077     using std::complex;
0078     using std::numeric_limits;
0079     constexpr const Real step = (numeric_limits<Real>::epsilon)();
0080     constexpr const Real inv_step = 1/(numeric_limits<Real>::epsilon)();
0081     return f(complex<Real>(x, step)).imag()*inv_step;
0082 }
0083
0084 namespace detail {
0085
0086    template <unsigned>
0087    struct fd_tag {};
0088
0089    template<class F, class Real>
0090    Real finite_difference_derivative(const F f, Real x, Real* error, const fd_tag<1>&)
0091    {
0092       using std::sqrt;
0093       using std::pow;
0094       using std::abs;
0095       using std::numeric_limits;
0096
0097       const Real eps = (numeric_limits<Real>::epsilon)();
0098       // Error bound ~eps^1/2
0099       // Note that this estimate of h differs from the best estimate by a factor of sqrt((|f(x)| + |f(x+h)|)/|f''(x)|).
0100       // Since this factor is invariant under the scaling f -> kf, then we are somewhat justified in approximating it by 1.
0101       // This approximation will get better as we move to higher orders of accuracy.
0102       Real h = 2 * sqrt(eps);
0103       h = detail::make_xph_representable(x, h);
0104
0105       Real yh = f(x + h);
0106       Real y0 = f(x);
0107       Real diff = yh - y0;
0108       if (error)
0109       {
0110          Real ym = f(x - h);
0111          Real ypph = abs(yh - 2 * y0 + ym) / h;
0112          // h*|f''(x)|*0.5 + (|f(x+h)+|f(x)|)*eps/h
0113          *error = ypph / 2 + (abs(yh) + abs(y0))*eps / h;
0114       }
0115       return diff / h;
0116    }
0117
0118    template<class F, class Real>
0119    Real finite_difference_derivative(const F f, Real x, Real* error, const fd_tag<2>&)
0120    {
0121       using std::sqrt;
0122       using std::pow;
0123       using std::abs;
0124       using std::numeric_limits;
0125
0126       const Real eps = (numeric_limits<Real>::epsilon)();
0127       // Error bound ~eps^2/3
0128       // See the previous discussion to understand determination of h and the error bound.
0129       // Series[(f[x+h] - f[x-h])/(2*h), {h, 0, 4}]
0130       Real h = pow(3 * eps, static_cast<Real>(1) / static_cast<Real>(3));
0131       h = detail::make_xph_representable(x, h);
0132
0133       Real yh = f(x + h);
0134       Real ymh = f(x - h);
0135       Real diff = yh - ymh;
0136       if (error)
0137       {
0138          Real yth = f(x + 2 * h);
0139          Real ymth = f(x - 2 * h);
0140          *error = eps * (abs(yh) + abs(ymh)) / (2 * h) + abs((yth - ymth) / 2 - diff) / (6 * h);
0141       }
0142
0143       return diff / (2 * h);
0144    }
0145
0146    template<class F, class Real>
0147    Real finite_difference_derivative(const F f, Real x, Real* error, const fd_tag<4>&)
0148    {
0149       using std::sqrt;
0150       using std::pow;
0151       using std::abs;
0152       using std::numeric_limits;
0153
0154       const Real eps = (numeric_limits<Real>::epsilon)();
0155       // Error bound ~eps^4/5
0156       Real h = pow(Real(11.25)*eps, static_cast<Real>(1) / static_cast<Real>(5));
0157       h = detail::make_xph_representable(x, h);
0158       Real ymth = f(x - 2 * h);
0159       Real yth = f(x + 2 * h);
0160       Real yh = f(x + h);
0161       Real ymh = f(x - h);
0162       Real y2 = ymth - yth;
0163       Real y1 = yh - ymh;
0164       if (error)
0165       {
0166          // Mathematica code to extract the remainder:
0167          // Series[(f[x-2*h]+ 8*f[x+h] - 8*f[x-h] - f[x+2*h])/(12*h), {h, 0, 7}]
0168          Real y_three_h = f(x + 3 * h);
0169          Real y_m_three_h = f(x - 3 * h);
0170          // Error from fifth derivative:
0171          *error = abs((y_three_h - y_m_three_h) / 2 + 2 * (ymth - yth) + 5 * (yh - ymh) / 2) / (30 * h);
0172          // Error from function evaluation:
0173          *error += eps * (abs(yth) + abs(ymth) + 8 * (abs(ymh) + abs(yh))) / (12 * h);
0174       }
0175       return (y2 + 8 * y1) / (12 * h);
0176    }
0177
0178    template<class F, class Real>
0179    Real finite_difference_derivative(const F f, Real x, Real* error, const fd_tag<6>&)
0180    {
0181       using std::sqrt;
0182       using std::pow;
0183       using std::abs;
0184       using std::numeric_limits;
0185
0186       const Real eps = (numeric_limits<Real>::epsilon)();
0187       // Error bound ~eps^6/7
0188       // Error: h^6f^(7)(x)/140 + 5|f(x)|eps/h
0189       Real h = pow(eps / 168, static_cast<Real>(1) / static_cast<Real>(7));
0190       h = detail::make_xph_representable(x, h);
0191
0192       Real yh = f(x + h);
0193       Real ymh = f(x - h);
0194       Real y1 = yh - ymh;
0195       Real y2 = f(x - 2 * h) - f(x + 2 * h);
0196       Real y3 = f(x + 3 * h) - f(x - 3 * h);
0197
0198       if (error)
0199       {
0200          // Mathematica code to generate fd scheme for 7th derivative:
0201          // Sum[(-1)^i*Binomial[7, i]*(f[x+(3-i)*h] + f[x+(4-i)*h])/2, {i, 0, 7}]
0202          // Mathematica to demonstrate that this is a finite difference formula for 7th derivative:
0203          // 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}]
0204          Real y7 = (f(x + 4 * h) - f(x - 4 * h) - 6 * y3 - 14 * y1 - 14 * y2) / 2;
0205          *error = abs(y7) / (140 * h) + 5 * (abs(yh) + abs(ymh))*eps / h;
0206       }
0207       return (y3 + 9 * y2 + 45 * y1) / (60 * h);
0208    }
0209
0210    template<class F, class Real>
0211    Real finite_difference_derivative(const F f, Real x, Real* error, const fd_tag<8>&)
0212    {
0213       using std::sqrt;
0214       using std::pow;
0215       using std::abs;
0216       using std::numeric_limits;
0217
0218       const Real eps = (numeric_limits<Real>::epsilon)();
0219       // Error bound ~eps^8/9.
0220       // In double precision, we only expect to lose two digits of precision while using this formula, at the cost of 8 function evaluations.
0221       // Error: h^8|f^(9)(x)|/630 + 7|f(x)|eps/h assuming 7 unstabilized additions.
0222       // Mathematica code to get the error:
0223       // 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}]
0224       // If we used Kahan summation, we could get the max error down to h^8|f^(9)(x)|/630 + |f(x)|eps/h.
0225       Real h = pow(Real(551.25)*eps, static_cast<Real>(1) / static_cast<Real>(9));
0226       h = detail::make_xph_representable(x, h);
0227
0228       Real yh = f(x + h);
0229       Real ymh = f(x - h);
0230       Real y1 = yh - ymh;
0231       Real y2 = f(x - 2 * h) - f(x + 2 * h);
0232       Real y3 = f(x + 3 * h) - f(x - 3 * h);
0233       Real y4 = f(x - 4 * h) - f(x + 4 * h);
0234
0235       Real tmp1 = 3 * y4 / 8 + 4 * y3;
0236       Real tmp2 = 21 * y2 + 84 * y1;
0237
0238       if (error)
0239       {
0240          // Mathematica code to generate fd scheme for 7th derivative:
0241          // Sum[(-1)^i*Binomial[9, i]*(f[x+(4-i)*h] + f[x+(5-i)*h])/2, {i, 0, 9}]
0242          // Mathematica to demonstrate that this is a finite difference formula for 7th derivative:
0243          // 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}]
0244          Real f9 = (f(x + 5 * h) - f(x - 5 * h)) / 2 + 4 * y4 + 27 * y3 / 2 + 24 * y2 + 21 * y1;
0245          *error = abs(f9) / (630 * h) + 7 * (abs(yh) + abs(ymh))*eps / h;
0246       }
0247       return (tmp1 + tmp2) / (105 * h);
0248    }
0249
0250    template<class F, class Real, class tag>
0251    Real finite_difference_derivative(const F, Real, Real*, const tag&)
0252    {
0253       // Always fails, but condition is template-arg-dependent so only evaluated if we get instantiated.
0254       static_assert(sizeof(Real) == 0, "Finite difference not implemented for this order: try 1, 2, 4, 6 or 8");
0255    }
0256
0257 }
0258
0259 template<class F, class Real, size_t order=6>
0260 inline Real finite_difference_derivative(const F f, Real x, Real* error = nullptr)
0261 {
0262    return detail::finite_difference_derivative(f, x, error, detail::fd_tag<order>());
0263 }
0264
0265 }}}  // namespaces
0266 #endif