3
\$\begingroup\$

Numerical differentiation is known to be ill-conditioned unless using a Chebyshev series, but this requires global information about the function and a priori knowledge of a compact domain on which the function will be evaluated. For this reason, simple finite differences are often useful. The unit roundoff gives a natural choice of step-see here and here for more details. I felt that this technique is better than the technique used in the GSL which requires specifying the stepsize, so I felt it should be given life in Boost.Math, esp. since spectral element methods are getting so popular. Regrettably, the pull-request was panned because the implementation is awkward, requiring a bunch of if checks on the order of accuracy and quite a bit of code duplication. Can it be cleaned up?

#ifndef BOOST_MATH_TOOLS_NUMERICAL_DIFFERENTIATION_HPP
#define BOOST_MATH_TOOLS_NUMERICAL_DIFFERENTIATION_HPP
#include <boost/math/constants/constants.hpp>

namespace boost { namespace math { namespace tools {
template<class F, class Real, size_t order=6>
Real finite_difference_derivative(const F f, Real x, Real* error = nullptr)
{
    static_assert(order == 1 || order == 2 || order == 4 || order == 6 || order == 8,
                  "Order of accuracy must be one of 1, 2, 4, 6, or 8.\n");
    // These using declarations allow for ADL so the boost.multiprecision types work.
    using std::sqrt;
    using std::pow;
    using std::abs;
    using std::max;
    using std::nextafter;
    using boost::math::constants::half;

    const Real eps = std::numeric_limits<Real>::epsilon();
    // TODO: static if in C++17!
    // Error bound ~eps^1/2
    if (order == 1)
    {
        // Note that this estimate of h differs from the best estimate by a factor of sqrt((|f(x)| + |f(x+h)|)/|f''(x)|).
        // Since this factor is invariant under the scaling f -> kf, then we are somewhat justified in approximating it by 1.
        // This approximation will get better as we move to higher orders of accuracy.
        Real h = 2*sqrt(eps);

        // Redefine h so that x + h is representable. Not using this trick leads to large error.
        // The compiler flag -ffast-math evaporates these operations . . .
        Real temp = x + h;
        h = temp - x;
        // Handle the case x + h == x:
        if (h == 0)
        {
            h = nextafter(x, std::numeric_limits<Real>::max()) - x;
        }
        Real yh = f(x+h);
        Real y0 = f(x);
        Real diff = yh - y0;
        if (error)
        {
            Real ym = f(x-h);
            Real ypph = abs(yh - 2*y0 + ym)/h;
            // h*|f''(x)|*0.5 + (|f(x+h)+|f(x)|)*eps/h, 
            *error = ypph*half<Real>() + (abs(yh) + abs(y0))*eps/h;
        }
        return diff/h;
    }

    // Error bound ~eps^2/3
    if (order == 2)
    {
        // See the previous discussion to understand determination of h and the error bound.
        // Series[(f[x+h] - f[x-h])/(2*h), {h, 0, 4}]
        Real h = pow(3*eps, boost::math::constants::third<Real>());
        Real temp = x + h;
        h = temp - x;
        if (h == 0)
        {
            h = nextafter(x, std::numeric_limits<Real>::max()) - x;
        }
        Real yh = f(x+h);
        Real ymh = f(x-h);
        Real diff = yh - ymh;
        if (error)
        {
            Real yth = f(x+2*h);
            Real ymth = f(x-2*h);
            *error = eps*(abs(yh) + abs(ymh))/(2*h) + abs((yth - ymth)*half<Real>() - diff)/(6*h);
        }

        return diff/(2*h);
    }

    // Error bound ~eps^4/5
    if (order == 4)
    {
        Real h = pow(11.25*eps, (Real) 1/ (Real) 5);
        Real temp = x + h;
        h = temp - x;
        if (h == 0)
        {
            h = nextafter(x, std::numeric_limits<Real>::max()) - x;
        }
        Real ymth = f(x-2*h);
        Real yth = f(x+2*h);
        Real yh = f(x+h);
        Real ymh = f(x-h);
        Real y2 = ymth - yth;
        Real y1 = yh - ymh;
        if (error)
        {
            // Mathematica code to extrace the remainder:
            // Series[(f[x-2*h]+ 8*f[x+h] - 8*f[x-h] - f[x+2*h])/(12*h), {h, 0, 7}]
            Real y_three_h = f(x+3*h);
            Real y_m_three_h = f(x-3*h);
            // Error from fifth derivative:
            *error = abs(half<Real>()*(y_three_h - y_m_three_h) + 2*(ymth - yth) + 5*half<Real>()*(yh - ymh) )/(30*h);
            // Error from function evaluation
            *error += eps*(abs(yth) + abs(ymth) + 8*(abs(ymh) + abs(yh)))/(12*h);
        }
        return (y2+8*y1)/(12*h);
    }

    // Error bound ~eps^6/7
    if (order == 6)
    {
        // Error: h^6f^(7)(x)/140 + 5|f(x)|eps/h
        Real h = pow(eps/168, (Real) 1/ (Real) 7);
        Real temp = x + h;
        h = temp - x;
        if (h == 0)
        {
            h = nextafter(x, std::numeric_limits<Real>::max()) - x;
        }

        Real yh = f(x+h);
        Real ymh = f(x-h);
        Real y1 = yh - ymh;
        Real y2 = f(x-2*h) - f(x + 2*h);
        Real y3 = f(x+3*h) - f(x - 3*h);

        if (error)
        {
            // Mathematica code to generate fd scheme for 7th derivative:
            // Sum[(-1)^i*Binomial[7, i]*(f[x+(3-i)*h] + f[x+(4-i)*h])/2, {i, 0, 7}]
            // Mathematica to demonstrate that this is a finite difference formula for 7th derivative:
            // 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}]
            Real y7 = half<Real>()*(f(x+4*h) - f(x-4*h) - 6*y3 - 14*y1 - 14*y2);
            *error = abs(y7)/(140*h) + 5*(abs(yh) + abs(ymh))*eps/h;
        }
        return (y3 + 9*y2 + 45*y1)/(60*h);
    }

    // Error bound ~eps^8/9.
    // In double precision, we only expect to lose two digits of precision while using this formula, at the cost of 8 function evaluations.
    if (order == 8)
    {
        // Error: h^8|f^(9)(x)|/630 + 7|f(x)|eps/h assuming 7 unstabilized additions.
        // 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, (Real)1 / (Real) 9);
        Real temp = x + h;
        h = temp - x;
        if (h == 0)
        {
            h = nextafter(x, std::numeric_limits<Real>::max()) - x;
        }

        Real yh = f(x+h);
        Real ymh = f(x-h);
        Real y1 = yh - ymh;
        Real y2 = f(x-2*h) - f(x + 2*h);
        Real y3 = f(x+3*h) - f(x - 3*h);
        Real y4 = f(x-4*h) - f(x + 4*h);

        Real tmp1 = 3*y4/8 + 4*y3;
        Real tmp2 = 21*y2 + 84*y1;

        if (error)
        {
            // Mathematica code to generate fd scheme for 7th derivative:
            // Sum[(-1)^i*Binomial[9, i]*(f[x+(4-i)*h] + f[x+(5-i)*h])/2, {i, 0, 9}]
            // Mathematica to demonstrate that this is a finite difference formula for 7th derivative:
            // 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}]
            Real f9 = (f(x+5*h) - f(x-5*h))*half<Real>() + 4*y4 + 27*half<Real>()*y3 + 24*y2 + 21*y1;
            *error = abs(f9)/(630*h) + 7*(abs(yh)+abs(ymh))*eps/h;
        }
        return (tmp1 + tmp2)/(105*h);
    }
}

}}}
#endif
\$\endgroup\$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Browse other questions tagged or ask your own question.