IT++ Logo
integration.h
Go to the documentation of this file.
1 
29 #ifndef INTEGRATION_H
30 #define INTEGRATION_H
31 
32 #include <limits>
35 #include <itpp/base/matfunc.h>
36 #include <itpp/base/specmat.h>
37 #include <itpp/itexports.h>
38 
39 namespace itpp
40 {
41 namespace details
42 {
44 template<typename Ftn>
45 double quadstep(Ftn f, double a, double b,
46  double fa, double fm, double fb, double is)
47 {
48  double Q, m, h, fml, fmr, i1, i2;
49 
50  m = (a + b) / 2;
51  h = (b - a) / 4;
52  fml = f(a + h);
53  fmr = f(b - h);
54 
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;
58 
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");
62  }
63  Q = i1;
64  return Q;
65  }
66  else {
67  Q = quadstep(f, a, m, fa, fml, fm, is) + quadstep(f, m, b, fm, fmr, fb, is);
68  }
69  return Q;
70 }
71 
73 template<typename Ftn>
74 double quadlstep(Ftn f, double a, double b,
75  double fa, double fb, double is)
76 {
77  static const double alpha = std::sqrt(2.0 / 3);
78  static const double beta = 1.0 / std::sqrt(5.0);
79 
80  double Q, h, m, mll, ml, mr, mrr, fmll, fml, fm, fmr, fmrr,
81  i1, i2;
82 
83  h = (b - a) / 2;
84  m = (a + b) / 2;
85 
86  mll = m - alpha * h;
87  ml = m - beta * h;
88  mr = m + beta * h;
89  mrr = m + alpha * h;
90 
91  fmll = f(mll);
92  fml = f(ml);
93  fm = f(m);
94  fmr = f(mr);
95  fmrr = f(mrr);
96 
97  i2 = (h / 6) * (fa + fb + 5 * (fml + fmr));
98  i1 = (h / 1470) * (77 * (fa + fb) + 432 * (fmll + fmrr) + 625 * (fml + fmr) + 672 * fm);
99 
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");
103  }
104  Q = i1;
105  return Q;
106  }
107  else {
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);
110  }
111  return Q;
112 }
113 }
114 
115 
122 
161 template <typename Ftn>
162 double quad(Ftn f, double a, double b,
163  double tol = std::numeric_limits<double>::epsilon())
164 {
165  vec x(3), y(3), yy(5);
166  double Q, fa, fm, fb, is;
167 
168  x = vec_3(a, (a + b) / 2, b);
169  y = apply_functor<double, Ftn>(f, x);
170  fa = y(0);
171  fm = y(1);
172  fb = y(2);
173  yy = apply_functor<double, Ftn>(f, a + vec(".9501 .2311 .6068 .4860 .8913")
174  * (b - a));
175  is = (b - a) / 8 * (sum(y) + sum(yy));
176 
177  if(is == 0.0)
178  is = b - a;
179 
180  is = is * tol / std::numeric_limits<double>::epsilon();
181  Q = details::quadstep(f, a, b, fa, fm, fb, is);
182 
183  return Q;
184 }
185 
222 ITPP_EXPORT double quad(double(*f)(double), double a, double b,
223  double tol = std::numeric_limits<double>::epsilon());
224 
263 template<typename Ftn>
264 double quadl(Ftn f, double a, double b,
265  double tol = std::numeric_limits<double>::epsilon())
266 {
267  static const double alpha = std::sqrt(2.0 / 3);
268  static const double beta = 1.0 / std::sqrt(5.0);
269 
270  double Q, m, h, x1, x2, x3, fa, fb, i1, i2, is, s, erri1, erri2, R;
271  vec x(13), y(13);
272  double tol2 = tol;
273 
274  m = (a + b) / 2;
275  h = (b - a) / 2;
276 
277  x1 = .942882415695480;
278  x2 = .641853342345781;
279  x3 = .236383199662150;
280  x(0) = a;
281  x(1) = m - x1 * h;
282  x(2) = m - alpha * h;
283  x(3) = m - x2 * h;
284  x(4) = m - beta * h;
285  x(5) = m - x3 * h;
286  x(6) = m;
287  x(7) = m + x3 * h;
288  x(8) = m + beta * h;
289  x(9) = m + x2 * h;
290  x(10) = m + alpha * h;
291  x(11) = m + x1 * h;
292  x(12) = b;
293 
294  y = apply_functor<double, Ftn>(f, x);
295 
296  fa = y(0);
297  fb = y(12);
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));
300 
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));
303 
304  s = sign(is);
305  if(s == 0.0)
306  s = 1;
307 
308  erri1 = std::abs(i1 - is);
309  erri2 = std::abs(i2 - is);
310 
311  R = 1;
312  if(erri2 != 0.0)
313  R = erri1 / erri2;
314 
315  if(R > 0 && R < 1)
316  tol2 = tol2 / R;
317 
318  is = s * std::abs(is) * tol2 / std::numeric_limits<double>::epsilon();
319  if(is == 0.0)
320  is = b - a;
321 
322  Q = details::quadlstep(f, a, b, fa, fb, is);
323 
324  return Q;
325 }
326 
363 ITPP_EXPORT double quadl(double(*f)(double), double a, double b,
364  double tol = std::numeric_limits<double>::epsilon());
365 
366 
368 
369 } // namespace itpp
370 
371 #endif // #ifndef INTEGRATION_H
Elementary mathematical functions - header file.
#define it_warning(s)
Display a warning message.
Definition: itassert.h:173
T sum(const Vec< T > &v)
Sum of all elements in the vector.
Definition: matfunc.h:59
double sign(double x)
Signum function.
Definition: elem_math.h:81
vec sqrt(const vec &x)
Square root of the elements.
Definition: elem_math.h:123
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.
Definition: integration.h:74
double quadstep(Ftn f, double a, double b, double fa, double fm, double fb, double is)
Simpson quadrature integration recursion step.
Definition: integration.h:45
Various functions on vectors and matrices - header file.
itpp namespace
Definition: itmex.h:37
double quadl(Ftn f, double a, double b, double tol=std::numeric_limits< double >::epsilon())
Definition: integration.h:264
bin abs(const bin &inbin)
absolute value of bin
Definition: binary.h:174
double quad(Ftn f, double a, double b, double tol=std::numeric_limits< double >::epsilon())
Definition: integration.h:162
Definitions of special vectors and matrices.

Generated on Tue Aug 17 2021 10:59:15 for IT++ by Doxygen 1.9.1