40 double quadstep(
double(*f)(
double),
double a,
double b,
41 double fa,
double fm,
double fb,
double is)
43 double Q, m, h, fml, fmr, i1, i2;
48 x = vec_2(a + h, b - h);
49 y = apply_function<double>(f, x);
53 i1 = h / 1.5 * (fa + 4 * fm + fb);
54 i2 = h / 3 * (fa + 4 * (fml + fmr) + 2 * fm + fb);
55 i1 = (16 * i2 - i1) / 15;
57 if ((is + (i1 - i2) == is) || (m <= a) || (b <= m)) {
58 if ((m <= a) || (b <= m)) {
59 it_warning(
"Interval contains no more machine number. Required tolerance may not be met");
65 Q =
quadstep(f, a, m, fa, fml, fm, is) +
quadstep(f, m, b, fm, fmr, fb, is);
71 double quad(
double(*f)(
double),
double a,
double b,
double tol)
73 vec x(3), y(3), yy(5);
74 double Q, fa, fm, fb, is;
76 x = vec_3(a, (a + b) / 2, b);
77 y = apply_function<double>(f, x);
81 yy = apply_function<double>(f, a + vec(
".9501 .2311 .6068 .4860 .8913")
83 is = (b - a) / 8 * (
sum(y) +
sum(yy));
88 is = is * tol / std::numeric_limits<double>::epsilon();
89 Q =
quadstep(f, a, b, fa, fm, fb, is);
98 double quadlstep(
double(*f)(
double),
double a,
double b,
99 double fa,
double fb,
double is)
101 double Q, h, m, alpha, beta, mll, ml, mr, mrr, fmll, fml, fm, fmr, fmrr,
119 y = apply_function<double>(f, x);
127 i2 = (h / 6) * (fa + fb + 5 * (fml + fmr));
128 i1 = (h / 1470) * (77 * (fa + fb) + 432 * (fmll + fmrr) + 625 * (fml + fmr) + 672 * fm);
130 if ((is + (i1 - i2) == is) || (mll <= a) || (b <= mrr)) {
131 if ((m <= a) || (b <= m)) {
132 it_warning(
"Interval contains no more machine number. Required tolerance may not be met");
138 Q =
quadlstep(f, a, mll, fa, fmll, is) +
quadlstep(f, mll, ml, fmll, fml, is) +
quadlstep(f, ml, m, fml, fm, is) +
139 quadlstep(f, m, mr, fm, fmr, is) +
quadlstep(f, mr, mrr, fmr, fmrr, is) +
quadlstep(f, mrr, b, fmrr, fb, is);
144 double quadl(
double(*f)(
double),
double a,
double b,
double tol)
146 double Q, m, h, alpha, beta, x1, x2, x3, fa, fb, i1, i2, is, s, erri1, erri2, R;
156 x1 = .942882415695480;
157 x2 = .641853342345781;
158 x3 = .236383199662150;
161 x(2) = m - alpha * h;
169 x(10) = m + alpha * h;
173 y = apply_function<double>(f, x);
177 i2 = (h / 6) * (y(0) + y(12) + 5 * (y(4) + y(8)));
178 i1 = (h / 1470) * (77 * (y(0) + y(12)) + 432 * (y(2) + y(10)) + 625 * (y(4) + y(8)) + 672 * y(6));
180 is = h * (.0158271919734802 * (y(0) + y(12)) + .0942738402188500 * (y(1) + y(11)) + .155071987336585 * (y(2) + y(10)) +
181 .188821573960182 * (y(3) + y(9)) + .199773405226859 * (y(4) + y(8)) + .224926465333340 * (y(5) + y(7)) + .242611071901408 * y(6));
197 is = s *
std::abs(is) * tol2 / std::numeric_limits<double>::epsilon();