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
Generated on Sat Jul 9 2011 15:21:30 for IT++ by Doxygen 1.7.4