src/nr/romberg.h

Go to the documentation of this file.
00001 //
00002 // C++ Interface: romberg
00003 //
00004 // Description: 
00005 //
00006 //
00007 // Author: Roman Geus <geus@maxwell>, (C) 2004
00008 //
00009 // Copyright: See COPYING file that comes with this distribution
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         // it = pow(2, n-2), n > 1
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;           // min. number of steps
00082     const int JMAX=20;          // max. number of steps
00083     const double EPS_REL=tol;   // relative accuracy desired
00084     const double EPS_ABS=1.0e-11;   // absolute accuracy desired
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;              // max. number of steps
00114     const int K=5;                  // 2*K is the order of the method
00115     const double EPS_REL=tol;       // relative accuracy desired
00116     const double EPS_ABS=1.0e-11;   // absolute accuracy desired
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 } // namespace NR
00146 
00147 #endif

Generated on Fri Oct 26 13:35:12 2007 for FEMAXX (Finite Element Maxwell Eigensolver) by  doxygen 1.4.7