00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013 #ifndef romberg_H
00014 #define romberg_H
00015
00016 #include <cassert>
00017 #include <cmath>
00018 #include <vector>
00019
00020 #include "rlog/rlog.h"
00021
00022 namespace NR {
00023
00024 using namespace std;
00025
00037 void polint(const vector<double>& xa, const vector<double>& ya, const double x, double &y, double &dy);
00038
00049 template <typename FunctionType>
00050 double trapzd(FunctionType& func, const double a, const double b, const int n) {
00051 double x, tnm, sum, del;
00052 static double s;
00053 int it, j;
00054
00055 if (n == 1) {
00056 return (s = 0.5*(b-a) * (func(a) + func(b)));
00057 } else {
00058
00059 for (it = 1, j = 1; j < n-1; ++ j)
00060 it <<= 1;
00061 tnm = it;
00062 del = (b-a)/tnm;
00063 x = a + 0.5*del;
00064 for (sum = 0.0, j = 0; j < it; ++ j, x += del)
00065 sum += func(x);
00066 s = 0.5*(s + (b-a)*sum/tnm);
00067 return s;
00068 }
00069 }
00070
00079 template <typename FunctionType>
00080 double qtrap(FunctionType& func, double a, double b, double tol, int& info) {
00081 const int JMIN=5;
00082 const int JMAX=20;
00083 const double EPS_REL=tol;
00084 const double EPS_ABS=1.0e-11;
00085 double s(0.0), olds;
00086
00087 info = 0;
00088 rAssert(JMIN > 0);
00089 for (int j = 0; j < JMAX; ++ j) {
00090 olds = s;
00091 s = trapzd<FunctionType>(func, a, b, j+1);
00092 rDebug("qtrap: j=%d, s=%e, ds=%e", j, s, std::fabs(s-olds));
00093 if (j > JMIN)
00094 if (std::fabs(s-olds) <= EPS_REL*fabs(olds) || std::fabs(s-olds) <= EPS_ABS)
00095 return s;
00096 }
00097 rWarning("qtrap not converged after %d steps with s=%e, ds=%e", JMAX, s, std::fabs(s-olds));
00098 info = 1;
00099 return s;
00100 }
00101
00111 template <typename FunctionType>
00112 double qromb(FunctionType& func, double a, double b, double tol, int& info) {
00113 const int JMAX=20;
00114 const int K=5;
00115 const double EPS_REL=tol;
00116 const double EPS_ABS=1.0e-11;
00117 const int JMAXP=JMAX+1;
00118
00119 double ss, dss;
00120 vector<double> s(JMAX), h(JMAXP);
00121 vector<double> s_t(K), h_t(K);
00122
00123 info = 0;
00124 h[0] = 1.0;
00125 for (int j = 1; j <= JMAX; ++ j) {
00126 s[j-1] = trapzd<FunctionType>(func, a, b, j);
00127 rDebug("qromb: j=%d, s=%e", j-1, s[j-1]);
00128 if (j >= K) {
00129 for (int i = 0; i < K; ++ i) {
00130 h_t[i] = h[j-K+i];
00131 s_t[i] = s[j-K+i];
00132 }
00133 polint(h_t, s_t, 0.0, ss, dss);
00134 rDebug("qromb: j=%d, ss=%e, dss=%e", j-1, ss, dss);
00135 if (std::fabs(dss) <= EPS_REL*fabs(ss) || std::fabs(dss) <= EPS_ABS)
00136 return ss;
00137 }
00138 h[j] = 0.25*h[j-1];
00139 }
00140 rWarning("qromb not converged after %d steps with ss=%e, dss=%e", JMAX, ss, dss);
00141 info = 1;
00142 return ss;
00143 }
00144
00145 }
00146
00147 #endif