NumericalDerivative.h

//
//  Copyright (C) 2017, Alexander Stukowski and Paul Erhart
//
//  This file is part of atomicrex.
//
//  Atomicrex is free software; you can redistribute it and/or modify
//  it under the terms of the GNU General Public License as published by
//  the Free Software Foundation; either version 3 of the License, or
//  (at your option) any later version.
//
//  Atomicrex is distributed in the hope that it will be useful,
//  but WITHOUT ANY WARRANTY; without even the implied warranty of
//  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
//  GNU General Public License for more details.
//
//  You should have received a copy of the GNU General Public License
//  along with this program.  If not, see <http://www.gnu.org/licenses/>.
//

#pragma once

#include "../Atomicrex.h"

namespace atomicrex {

class NumericalDerivative
{
public:

    double compute(double x, double& error, double epsilon = 1e-5) {

        double round, trunc;
        double r_0 = centralDerivative(x, epsilon, round, trunc);
        error = round + trunc;
        if(round < trunc && (round > 0 && trunc > 0)) {
            // Compute an optimized step size to minimize the total error,
            // using the scaling of the truncation error (O(h^2)) and
            // rounding error (O(1/h)).
            double h_opt = epsilon * pow(round / (2.0 * trunc), 1.0 / 3.0);
            double round_opt, trunc_opt;
            double r_opt = centralDerivative(x, h_opt, round_opt, trunc_opt);
            double error_opt = round_opt + trunc_opt;

            // Check that the new error is smaller, and that the new derivative
            // is consistent with the error bounds of the original estimate.
            if(error_opt < error && fabs(r_opt - r_0) < 4.0 * error) {
                r_0 = r_opt;
                error = error_opt;
            }
        }

        return r_0;
    }

protected:

    virtual double evaluate(double x) = 0;

    double centralDerivative(double x, double h, double& abserr_trunc, double& abserr_round) {

        // Compute the error using the difference between the 5-point and
        // the 3-point rule (-h,0,+h). Again the central point is not used.

        double fm1 = evaluate(x - h);
        double fp1 = evaluate(x + h);

        double fmh = evaluate(x - h/2);
        double fph = evaluate(x + h/2);

        double r3 = 0.5 * (fp1 - fm1);
        double r5 = (4.0 / 3.0) * (fph - fmh) - (1.0 / 3.0) * r3;

        const double CENTRAL_DBL_EPSILON = 2.2204460492503131e-16;
        double e3 = (fabs(fp1) + fabs(fm1)) * CENTRAL_DBL_EPSILON;
        double e5 = 2.0 * (fabs(fph) + fabs(fmh)) * CENTRAL_DBL_EPSILON + e3;

        // The truncation error in the r5 approximation itself is O(h^4).
        // However, for safety, we estimate the error from r5-r3, which is
        // O(h^2).  By scaling h we will minimize this estimated error, not
        // the actual truncation error in r5.

        double result = r5 / h;
        abserr_trunc = fabs((r5 - r3) / h); // Estimated truncation error O(h^2)
        abserr_round = fabs(e5 / h);        // Rounding error (cancellations)

        return result;
    }
};

} // End of namespace