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
Generated on Sat Apr 19 10:57:51 2008 for IT++ by Doxygen 1.5.5