src/nr/romberg.h

Go to the documentation of this file.

//
// C++ Interface: romberg
//
// Description: 
//
//
// Author: Roman Geus <geus@maxwell>, (C) 2004
//
// Copyright: See COPYING file that comes with this distribution
//
//

#ifndef romberg_H
#define romberg_H

#include <cassert>
#include <cmath>
#include <vector>

#include "rlog/rlog.h"

namespace NR {

using namespace std;

void polint(const vector<double>& xa, const vector<double>& ya, const double x, double &y, double &dy);

template <typename FunctionType>
double trapzd(FunctionType& func, const double a, const double b, const int n) {
    double x, tnm, sum, del;
    static double s;
    int it, j;

    if (n == 1) {
        return (s = 0.5*(b-a) * (func(a) + func(b)));
    } else {
        // it = pow(2, n-2), n > 1
        for (it = 1, j = 1; j < n-1; ++ j)
            it <<= 1;
        tnm = it;
        del = (b-a)/tnm;
        x = a + 0.5*del;
        for (sum = 0.0, j = 0; j < it; ++ j, x += del)
            sum += func(x);
        s = 0.5*(s + (b-a)*sum/tnm);
        return s;
    }
}

template <typename FunctionType>
double qtrap(FunctionType& func, double a, double b, double tol, int& info) {
    const int JMIN=5;           // min. number of steps
    const int JMAX=20;          // max. number of steps
    const double EPS_REL=tol;   // relative accuracy desired
    const double EPS_ABS=1.0e-11;   // absolute accuracy desired
    double s(0.0), olds;
    
    info = 0;
    rAssert(JMIN > 0);
    for (int j = 0; j < JMAX; ++ j) {
        olds = s;
        s = trapzd<FunctionType>(func, a, b, j+1);
        rDebug("qtrap: j=%d, s=%e, ds=%e", j, s, std::fabs(s-olds));
        if (j > JMIN)
            if (std::fabs(s-olds) <= EPS_REL*fabs(olds) || std::fabs(s-olds) <= EPS_ABS)
                return s;
    }
    rWarning("qtrap not converged after %d steps with s=%e, ds=%e", JMAX, s, std::fabs(s-olds));
    info = 1;
    return s;
}

template <typename FunctionType>
double qromb(FunctionType& func, double a, double b, double tol, int& info) {
    const int JMAX=20;              // max. number of steps
    const int K=5;                  // 2*K is the order of the method
    const double EPS_REL=tol;       // relative accuracy desired
    const double EPS_ABS=1.0e-11;   // absolute accuracy desired
    const int JMAXP=JMAX+1;

    double ss, dss;
    vector<double> s(JMAX), h(JMAXP);
    vector<double> s_t(K), h_t(K);
    
    info = 0;
    h[0] = 1.0;
    for (int j = 1; j <= JMAX; ++ j) {
        s[j-1] = trapzd<FunctionType>(func, a, b, j);
        rDebug("qromb: j=%d, s=%e", j-1, s[j-1]);
        if (j >= K) {
            for (int i = 0; i < K; ++ i) {
                h_t[i] = h[j-K+i];
                s_t[i] = s[j-K+i];
            }
            polint(h_t, s_t, 0.0, ss, dss);
            rDebug("qromb: j=%d, ss=%e, dss=%e", j-1, ss, dss);
            if (std::fabs(dss) <= EPS_REL*fabs(ss) || std::fabs(dss) <= EPS_ABS)
                return ss;
        }
        h[j] = 0.25*h[j-1];
    }
    rWarning("qromb not converged after %d steps with ss=%e, dss=%e", JMAX, ss, dss);
    info = 1;
    return ss;
}

} // namespace NR

#endif

---