IT++ Logo
integration.cpp
Go to the documentation of this file.
00001 
00029 #include <itpp/base/math/integration.h>
00030 #include <itpp/base/math/elem_math.h>
00031 #include <itpp/base/help_functions.h>
00032 #include <itpp/base/matfunc.h>
00033 #include <itpp/base/specmat.h>
00034 
00035 
00036 namespace itpp
00037 {
00038 
00040 double quadstep(double(*f)(double), double a, double b,
00041                 double fa, double fm, double fb, double is)
00042 {
00043   double Q, m, h, fml, fmr, i1, i2;
00044   vec x(2), y(2);
00045 
00046   m = (a + b) / 2;
00047   h = (b - a) / 4;
00048   x = vec_2(a + h, b - h);
00049   y = apply_function<double>(f, x);
00050   fml = y(0);
00051   fmr = y(1);
00052 
00053   i1 = h / 1.5 * (fa + 4 * fm + fb);
00054   i2 = h / 3 * (fa + 4 * (fml + fmr) + 2 * fm + fb);
00055   i1 = (16 * i2 - i1) / 15;
00056 
00057   if ((is + (i1 - i2) == is) || (m <= a) || (b <= m)) {
00058     if ((m <= a) || (b <= m)) {
00059       it_warning("Interval contains no more machine number. Required tolerance may not be met");
00060     }
00061     Q = i1;
00062     return Q;
00063   }
00064   else {
00065     Q = quadstep(f, a, m, fa, fml, fm, is) + quadstep(f, m, b, fm, fmr, fb, is);
00066   }
00067   return Q;
00068 }
00069 
00070 
00071 double quad(double(*f)(double), double a, double b, double tol)
00072 {
00073   vec x(3), y(3), yy(5);
00074   double Q, fa, fm, fb, is;
00075 
00076   x = vec_3(a, (a + b) / 2, b);
00077   y = apply_function<double>(f, x);
00078   fa = y(0);
00079   fm = y(1);
00080   fb = y(2);
00081   yy = apply_function<double>(f, a + vec(".9501 .2311 .6068 .4860 .8913")
00082                               * (b - a));
00083   is = (b - a) / 8 * (sum(y) + sum(yy));
00084 
00085   if (is == 0.0)
00086     is = b - a;
00087 
00088   is = is * tol / std::numeric_limits<double>::epsilon();
00089   Q = quadstep(f, a, b, fa, fm, fb, is);
00090 
00091   return Q;
00092 }
00093 
00094 
00095 //--------------------- quadl() ----------------------------------------
00096 
00098 double quadlstep(double(*f)(double), double a, double b,
00099                  double fa, double fb, double is)
00100 {
00101   double Q, h, m, alpha, beta, mll, ml, mr, mrr, fmll, fml, fm, fmr, fmrr,
00102   i1, i2;
00103   vec x(5), y(5);
00104 
00105   h = (b - a) / 2;
00106   m = (a + b) / 2;
00107   alpha = std::sqrt(2.0 / 3);
00108   beta = 1.0 / std::sqrt(5.0);
00109   mll = m - alpha * h;
00110   ml = m - beta * h;
00111   mr = m + beta * h;
00112   mrr = m + alpha * h;
00113   x(0) = mll;
00114   x(1) = ml;
00115   x(2) = m;
00116   x(3) = mr;
00117   x(4) = mrr;
00118 
00119   y = apply_function<double>(f, x);
00120 
00121   fmll = y(0);
00122   fml = y(1);
00123   fm = y(2);
00124   fmr = y(3);
00125   fmrr = y(4);
00126 
00127   i2 = (h / 6) * (fa + fb + 5 * (fml + fmr));
00128   i1 = (h / 1470) * (77 * (fa + fb) + 432 * (fmll + fmrr) + 625 * (fml + fmr) + 672 * fm);
00129 
00130   if ((is + (i1 - i2) == is) || (mll <= a) || (b <= mrr)) {
00131     if ((m <= a) || (b <= m)) {
00132       it_warning("Interval contains no more machine number. Required tolerance may not be met");
00133     }
00134     Q = i1;
00135     return Q;
00136   }
00137   else {
00138     Q = quadlstep(f, a, mll, fa, fmll, is) + quadlstep(f, mll, ml, fmll, fml, is) + quadlstep(f, ml, m, fml, fm, is) +
00139         quadlstep(f, m, mr, fm, fmr, is) + quadlstep(f, mr, mrr, fmr, fmrr, is) + quadlstep(f, mrr, b, fmrr, fb, is);
00140   }
00141   return Q;
00142 }
00143 
00144 double quadl(double(*f)(double), double a, double b, double tol)
00145 {
00146   double Q, m, h, alpha, beta, x1, x2, x3, fa, fb, i1, i2, is, s, erri1, erri2, R;
00147   vec x(13), y(13);
00148   double tol2 = tol;
00149 
00150   m = (a + b) / 2;
00151   h = (b - a) / 2;
00152 
00153   alpha = std::sqrt(2.0 / 3);
00154   beta = 1.0 / std::sqrt(5.0);
00155 
00156   x1 = .942882415695480;
00157   x2 = .641853342345781;
00158   x3 = .236383199662150;
00159   x(0) = a;
00160   x(1) = m - x1 * h;
00161   x(2) = m - alpha * h;
00162   x(3) = m - x2 * h;
00163   x(4) = m - beta * h;
00164   x(5) = m - x3 * h;
00165   x(6) = m;
00166   x(7) = m + x3 * h;
00167   x(8) = m + beta * h;
00168   x(9) = m + x2 * h;
00169   x(10) = m + alpha * h;
00170   x(11) = m + x1 * h;
00171   x(12) = b;
00172 
00173   y = apply_function<double>(f, x);
00174 
00175   fa = y(0);
00176   fb = y(12);
00177   i2 = (h / 6) * (y(0) + y(12) + 5 * (y(4) + y(8)));
00178   i1 = (h / 1470) * (77 * (y(0) + y(12)) + 432 * (y(2) + y(10)) + 625 * (y(4) + y(8)) + 672 * y(6));
00179 
00180   is = h * (.0158271919734802 * (y(0) + y(12)) + .0942738402188500 * (y(1) + y(11)) + .155071987336585 * (y(2) + y(10)) +
00181             .188821573960182 * (y(3) + y(9)) + .199773405226859 * (y(4) + y(8)) + .224926465333340 * (y(5) + y(7)) + .242611071901408 * y(6));
00182 
00183   s = sign(is);
00184   if (s == 0.0)
00185     s = 1;
00186 
00187   erri1 = std::abs(i1 - is);
00188   erri2 = std::abs(i2 - is);
00189 
00190   R = 1;
00191   if (erri2 != 0.0)
00192     R = erri1 / erri2;
00193 
00194   if (R > 0 && R < 1)
00195     tol2 = tol2 / R;
00196 
00197   is = s * std::abs(is) * tol2 / std::numeric_limits<double>::epsilon();
00198   if (is == 0.0)
00199     is = b - a;
00200 
00201   Q = quadlstep(f, a, b, fa, fb, is);
00202 
00203   return Q;
00204 }
00205 
00206 
00207 } // namespace itpp
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines
SourceForge Logo

Generated on Wed Jul 27 2011 16:27:04 for IT++ by Doxygen 1.7.4