IT++ Logo

integration.cpp

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

Generated on Sat Apr 19 10:41:11 2008 for IT++ by Doxygen 1.5.5