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());
Elementary mathematical functions - header file.
#define it_warning(s)
Display a warning message.
T sum(const Vec< T > &v)
Sum of all elements in the vector.
double sign(double x)
Signum function.
vec sqrt(const vec &x)
Square root of the elements.
Help functions to make functions with vec and mat as arguments.
double quadlstep(Ftn f, double a, double b, double fa, double fb, double is)
Adaptive Lobatto quadrature integration recursion step.
double quadstep(Ftn f, double a, double b, double fa, double fm, double fb, double is)
Simpson quadrature integration recursion step.
Various functions on vectors and matrices - header file.
double quadl(Ftn f, double a, double b, double tol=std::numeric_limits< double >::epsilon())
bin abs(const bin &inbin)
absolute value of bin
double quad(Ftn f, double a, double b, double tol=std::numeric_limits< double >::epsilon())
Definitions of special vectors and matrices.