NumericalIntegration.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 NumericalIntegration
{
public:

    template<typename DifferentialForm>
    double integrate(double a, double b, DifferentialForm diffForm) {
        if(a > b)
            throw std::runtime_error(boost::str(boost::format("Left boundary larger than right boundary in NumericalIntegration::integrate(): %1% > %2%.") % a % b));
        return simpson(a, b, diffForm);
    };

    virtual double integrate(double a, double b) {
        return integrate(a, b, differentialFormDefault);
    }

    virtual double integrateSpherical(double a, double b) {
        return integrate(a, b, differentialFormSpherical);
    }

protected:

    virtual double evaluateInternal(double x) = 0;

private:

    static double differentialFormDefault(double r, double f) {
        return f;
    }

    static double differentialFormSpherical(double r, double f) {
        return square(r) * f;
    }

private:

    /**********************************************************************
     * Returns the integral of the function func from a to b. The parameters
     * epsilon can be set to the desired fractional accuracy and JMAX so that 2
     * to the power JMAX-1 is the maximum allowed number of steps.
     * Integration is performed by Simpson's rule.
     **********************************************************************/
    template<typename DifferentialForm>
    double simpson(double a, double b, DifferentialForm diffForm, double eps = FLOATTYPE_EPSILON, int jmax = 20)
    {
        double s, s2, st, ost = 0.0, os = 0.0;
        for(int j = 1; j <= jmax; j++) {
            st = trapezoid(a, b, j, s2, diffForm);
            s = (4.0 * st - ost) / 3.0;
            if(j > 5) // Avoid spurious early convergence.
                if(fabs(s - os) < eps*fabs(os) || (s == 0.0 && os == 0.0)) return s;
            os = s;
            ost = st;
        }
        throw std::runtime_error("Exceeded maximum number of iterations in routine NumericalIntegration::simpson().");
        return 0.0; // Never get here.
    };

    /**********************************************************************
     * This routine computes the nth stage of refinement of an extended
     * trapezoidal rule. func is input as a pointer to the function to
     * be integrated between limits a and b, also input. When called with
     * n=1, the routine returns the crudest estimate of a int_a^b f(x)dx.
     * Subsequent calls with n=2,3,... (in that sequential order) will
     * improve the accuracy by adding 2^(n-2) additional interior points.
     **********************************************************************/
    template<typename DifferentialForm>
    double trapezoid(double a, double b, int n, double& s, DifferentialForm diffForm)
    {
        if(n == 1) {
            return (s = 0.5*(b-a)*(diffForm(a, evaluateInternal(a)) + diffForm(b, evaluateInternal(b))));
        }
        else {
            int it = 1;
            for(int j = 1; j < n-1; j++) it <<= 1;
            double tnm = it;
            double del = (b-a)/tnm; // This is the spacing of the points to be added.
            double x = a + 0.5 * del;
            double sum = 0;
            for(int j = 1; j <= it; j++, x+=del) sum += diffForm(x, evaluateInternal(x));
            s = 0.5 * (s + (b-a) * sum / tnm); // This replaces s by its refined value.
            return s;
        }
    }
};

} // End of namespace