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