src/nr/adaptsim.h

Go to the documentation of this file.

//
// C++ Implementation: adaptsim
//
// Description: Numerically evaluate integral using adaptive Simpson rule
//
//
// Author: Roman Geus <geus@maxwell>, (C) 2004
//
// Copyright: See COPYING file that comes with this distribution
//
//

#include <limits>

#ifndef adaptsim_H
#define adaptsim_H

namespace NR {

template <typename FunctionType, typename ScalarType>
ScalarType adaptsim(FunctionType& func, ScalarType a, ScalarType b, ScalarType tol, int& info) {
    const ScalarType mach_eps = std::numeric_limits<ScalarType>::epsilon();

    info = 0;
    if (tol < mach_eps)
        tol = mach_eps;

    double sign = 1.0;
    if (a > b) {
        sign = -1.0;
        double t = a;
        a = b;
        b = t;
    }   

    ScalarType fa = func(a);
    ScalarType fm = func(0.5*(a+b));
    ScalarType fb = func(b);

    // Compute rough estimate of modulus of integral is
    ScalarType mc_pos[] = { 0.9501, 0.2311, 0.6068, 0.4860, 0.8913 };
    ScalarType sum_y = fa + fm + fb;
    for (int i = 0; i < 5; ++ i)
        sum_y += func(a + mc_pos[i] * (b-a));
    ScalarType is = (b - a)/8.0 * sum_y;
    if (is == 0.0)
        is = b-a;

    // rDebug("adaptsim: is=%e", is);
    is *= tol/mach_eps;
    return sign*adaptsimstp(func, a, b, fa, fm, fb, is, info);
}

template <typename FunctionType, typename ScalarType>
ScalarType adaptsimstp(FunctionType& func,
                       ScalarType a, ScalarType b,
                       ScalarType fa, ScalarType fm, ScalarType fb,
                       ScalarType is,
                       int& info) {
    const ScalarType eps = std::numeric_limits<ScalarType>::epsilon() * fabs(is);
    ScalarType Q;
    ScalarType m = (a + b)/2;
    ScalarType h = (b - a)/4;
    ScalarType fml = func(a + h);
    ScalarType fmr = func(b - h);

    // Simpson's rule for h: i1 is 3-point formula
    ScalarType i1 = h/1.5 * (fa + 4*fm + fb);
    // Simpson's rule for h/2
    ScalarType i2 = h/3 * (fa + 4*(fml + fmr) + 2*fm + fb);
    // Bode's rule: i1 is 5-point formula
    i1 = (16*i2 - i1)/15;

    if ((fabs(i1-i2) < eps) || (m <= a) || (b <= m)) {
        if (((m <= a) || (b <= m)) && (info == 0)) {
            rWarning("Interval contains no more machine number. Required tolerance may not be met.");
            info = 1;
        }
        Q = i1;
        // rDebug("adaptsim: a=%e b-a=%e Q=%e", a, b-a, Q);
    } else {
        Q = adaptsimstp(func, a, m, fa, fml, fm, is, info) +
            adaptsimstp(func, m, b, fm, fmr, fb, is, info);
    }
    return Q;
}

}

#endif

---