37 #include <itpp/itexports.h> 44 template<
typename Ftn>
46 double fa,
double fm,
double fb,
double is)
48 double Q, m, h, fml, fmr, i1, i2;
55 i1 = h / 1.5 * (fa + 4 * fm + fb);
56 i2 = h / 3 * (fa + 4 * (fml + fmr) + 2 * fm + fb);
57 i1 = (16 * i2 - i1) / 15;
59 if((is + (i1 - i2) == is) || (m <= a) || (b <= m)) {
60 if((m <= a) || (b <= m)) {
61 it_warning(
"Interval contains no more machine number. Required tolerance may not be met");
67 Q =
quadstep(f, a, m, fa, fml, fm, is) +
quadstep(f, m, b, fm, fmr, fb, is);
73 template<
typename Ftn>
75 double fa,
double fb,
double is)
77 static const double alpha =
std::sqrt(2.0 / 3);
78 static const double beta = 1.0 /
std::sqrt(5.0);
80 double Q, h, m, mll, ml, mr, mrr, fmll, fml, fm, fmr, fmrr,
97 i2 = (h / 6) * (fa + fb + 5 * (fml + fmr));
98 i1 = (h / 1470) * (77 * (fa + fb) + 432 * (fmll + fmrr) + 625 * (fml + fmr) + 672 * fm);
100 if((is + (i1 - i2) == is) || (mll <= a) || (b <= mrr)) {
101 if((m <= a) || (b <= m)) {
102 it_warning(
"Interval contains no more machine number. Required tolerance may not be met");
108 Q =
quadlstep(f, a, mll, fa, fmll, is) +
quadlstep(f, mll, ml, fmll, fml, is) +
quadlstep(f, ml, m, fml, fm, is) +
109 quadlstep(f, m, mr, fm, fmr, is) +
quadlstep(f, mr, mrr, fmr, fmrr, is) +
quadlstep(f, mrr, b, fmrr, fb, is);
161 template <
typename Ftn>
162 double quad(Ftn f,
double a,
double b,
163 double tol = std::numeric_limits<double>::epsilon())
165 vec x(3), y(3), yy(5);
166 double Q, fa, fm, fb, is;
168 x = vec_3(a, (a + b) / 2, b);
169 y = apply_functor<double, Ftn>(f, x);
173 yy = apply_functor<double, Ftn>(f, a + vec(
".9501 .2311 .6068 .4860 .8913")
175 is = (b - a) / 8 * (
sum(y) +
sum(yy));
180 is = is * tol / std::numeric_limits<double>::epsilon();
181 Q = details::quadstep(f, a, b, fa, fm, fb, is);
222 ITPP_EXPORT
double quad(
double(*f)(
double),
double a,
double b,
223 double tol = std::numeric_limits<double>::epsilon());
263 template<
typename Ftn>
264 double quadl(Ftn f,
double a,
double b,
265 double tol = std::numeric_limits<double>::epsilon())
267 static const double alpha =
std::sqrt(2.0 / 3);
268 static const double beta = 1.0 /
std::sqrt(5.0);
270 double Q, m, h, x1, x2, x3, fa, fb, i1, i2, is, s, erri1, erri2, R;
277 x1 = .942882415695480;
278 x2 = .641853342345781;
279 x3 = .236383199662150;
282 x(2) = m - alpha * h;
290 x(10) = m + alpha * h;
294 y = apply_functor<double, Ftn>(f, x);
298 i2 = (h / 6) * (y(0) + y(12) + 5 * (y(4) + y(8)));
299 i1 = (h / 1470) * (77 * (y(0) + y(12)) + 432 * (y(2) + y(10)) + 625 * (y(4) + y(8)) + 672 * y(6));
301 is = h * (.0158271919734802 * (y(0) + y(12)) + .0942738402188500 * (y(1) + y(11)) + .155071987336585 * (y(2) + y(10)) +
302 .188821573960182 * (y(3) + y(9)) + .199773405226859 * (y(4) + y(8)) + .224926465333340 * (y(5) + y(7)) + .242611071901408 * y(6));
318 is = s *
std::abs(is) * tol2 / std::numeric_limits<double>::epsilon();
322 Q = details::quadlstep(f, a, b, fa, fb, is);
363 ITPP_EXPORT
double quadl(
double(*f)(
double),
double a,
double b,
364 double tol = std::numeric_limits<double>::epsilon());
371 #endif // #ifndef INTEGRATION_H Various functions on vectors and matrices - header file.
double quadstep(Ftn f, double a, double b, double fa, double fm, double fb, double is)
Simpson quadrature integration recursion step.
double quadlstep(Ftn f, double a, double b, double fa, double fb, double is)
Adaptive Lobatto quadrature integration recursion step.
Help functions to make functions with vec and mat as arguments.
T sum(const Vec< T > &v)
Sum of all elements in the vector.
double quad(Ftn f, double a, double b, double tol=std::numeric_limits< double >::epsilon())
double quadl(Ftn f, double a, double b, double tol=std::numeric_limits< double >::epsilon())
Definitions of special vectors and matrices.
#define it_warning(s)
Display a warning message.
double sign(double x)
Signum function.
vec sqrt(const vec &x)
Square root of the elements.
bin abs(const bin &inbin)
absolute value of bin
Elementary mathematical functions - header file.