Go to the documentation of this file.
```00001 //
00003 //
00004 // Description: Numerically evaluate integral using adaptive Simpson rule
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 #include <limits>
00014
00017
00018 namespace NR {
00019
00047 template <typename FunctionType, typename ScalarType>
00048 ScalarType adaptsim(FunctionType& func, ScalarType a, ScalarType b, ScalarType tol, int& info) {
00049     const ScalarType mach_eps = std::numeric_limits<ScalarType>::epsilon();
00050
00051     info = 0;
00052     if (tol < mach_eps)
00053         tol = mach_eps;
00054
00055     double sign = 1.0;
00056     if (a > b) {
00057         sign = -1.0;
00058         double t = a;
00059         a = b;
00060         b = t;
00061     }
00062
00063     ScalarType fa = func(a);
00064     ScalarType fm = func(0.5*(a+b));
00065     ScalarType fb = func(b);
00066
00067     // Compute rough estimate of modulus of integral is
00068     ScalarType mc_pos[] = { 0.9501, 0.2311, 0.6068, 0.4860, 0.8913 };
00069     ScalarType sum_y = fa + fm + fb;
00070     for (int i = 0; i < 5; ++ i)
00071         sum_y += func(a + mc_pos[i] * (b-a));
00072     ScalarType is = (b - a)/8.0 * sum_y;
00073     if (is == 0.0)
00074         is = b-a;
00075
00077     is *= tol/mach_eps;
00078     return sign*adaptsimstp(func, a, b, fa, fm, fb, is, info);
00079 }
00080
00093 template <typename FunctionType, typename ScalarType>
00095                        ScalarType a, ScalarType b,
00096                        ScalarType fa, ScalarType fm, ScalarType fb,
00097                        ScalarType is,
00098                        int& info) {
00099     const ScalarType eps = std::numeric_limits<ScalarType>::epsilon() * fabs(is);
00100     ScalarType Q;
00101     ScalarType m = (a + b)/2;
00102     ScalarType h = (b - a)/4;
00103     ScalarType fml = func(a + h);
00104     ScalarType fmr = func(b - h);
00105
00106     // Simpson's rule for h: i1 is 3-point formula
00107     ScalarType i1 = h/1.5 * (fa + 4*fm + fb);
00108     // Simpson's rule for h/2
00109     ScalarType i2 = h/3 * (fa + 4*(fml + fmr) + 2*fm + fb);
00110     // Bode's rule: i1 is 5-point formula
00111     i1 = (16*i2 - i1)/15;
00112
00113     if ((fabs(i1-i2) < eps) || (m <= a) || (b <= m)) {
00114         if (((m <= a) || (b <= m)) && (info == 0)) {
00115             rWarning("Interval contains no more machine number. Required tolerance may not be met.");
00116             info = 1;
00117         }
00118         Q = i1;
00119         // rDebug("adaptsim: a=%e b-a=%e Q=%e", a, b-a, Q);
00120     } else {
00121         Q = adaptsimstp(func, a, m, fa, fml, fm, is, info) +
00122             adaptsimstp(func, m, b, fm, fmr, fb, is, info);
00123     }
00124     return Q;
00125 }
00126
00127 }
00128
00129 #endif
```

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