IT++ Logo
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
integration.cpp
Go to the documentation of this file.
1 
32 #include <itpp/base/matfunc.h>
33 #include <itpp/base/specmat.h>
34 
35 
36 namespace itpp
37 {
38 
40 double quadstep(double(*f)(double), double a, double b,
41  double fa, double fm, double fb, double is)
42 {
43  double Q, m, h, fml, fmr, i1, i2;
44  vec x(2), y(2);
45 
46  m = (a + b) / 2;
47  h = (b - a) / 4;
48  x = vec_2(a + h, b - h);
49  y = apply_function<double>(f, x);
50  fml = y(0);
51  fmr = y(1);
52 
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;
56 
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");
60  }
61  Q = i1;
62  return Q;
63  }
64  else {
65  Q = quadstep(f, a, m, fa, fml, fm, is) + quadstep(f, m, b, fm, fmr, fb, is);
66  }
67  return Q;
68 }
69 
70 
71 double quad(double(*f)(double), double a, double b, double tol)
72 {
73  vec x(3), y(3), yy(5);
74  double Q, fa, fm, fb, is;
75 
76  x = vec_3(a, (a + b) / 2, b);
77  y = apply_function<double>(f, x);
78  fa = y(0);
79  fm = y(1);
80  fb = y(2);
81  yy = apply_function<double>(f, a + vec(".9501 .2311 .6068 .4860 .8913")
82  * (b - a));
83  is = (b - a) / 8 * (sum(y) + sum(yy));
84 
85  if (is == 0.0)
86  is = b - a;
87 
88  is = is * tol / std::numeric_limits<double>::epsilon();
89  Q = quadstep(f, a, b, fa, fm, fb, is);
90 
91  return Q;
92 }
93 
94 
95 //--------------------- quadl() ----------------------------------------
96 
98 double quadlstep(double(*f)(double), double a, double b,
99  double fa, double fb, double is)
100 {
101  double Q, h, m, alpha, beta, mll, ml, mr, mrr, fmll, fml, fm, fmr, fmrr,
102  i1, i2;
103  vec x(5), y(5);
104 
105  h = (b - a) / 2;
106  m = (a + b) / 2;
107  alpha = std::sqrt(2.0 / 3);
108  beta = 1.0 / std::sqrt(5.0);
109  mll = m - alpha * h;
110  ml = m - beta * h;
111  mr = m + beta * h;
112  mrr = m + alpha * h;
113  x(0) = mll;
114  x(1) = ml;
115  x(2) = m;
116  x(3) = mr;
117  x(4) = mrr;
118 
119  y = apply_function<double>(f, x);
120 
121  fmll = y(0);
122  fml = y(1);
123  fm = y(2);
124  fmr = y(3);
125  fmrr = y(4);
126 
127  i2 = (h / 6) * (fa + fb + 5 * (fml + fmr));
128  i1 = (h / 1470) * (77 * (fa + fb) + 432 * (fmll + fmrr) + 625 * (fml + fmr) + 672 * fm);
129 
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");
133  }
134  Q = i1;
135  return Q;
136  }
137  else {
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);
140  }
141  return Q;
142 }
143 
144 double quadl(double(*f)(double), double a, double b, double tol)
145 {
146  double Q, m, h, alpha, beta, x1, x2, x3, fa, fb, i1, i2, is, s, erri1, erri2, R;
147  vec x(13), y(13);
148  double tol2 = tol;
149 
150  m = (a + b) / 2;
151  h = (b - a) / 2;
152 
153  alpha = std::sqrt(2.0 / 3);
154  beta = 1.0 / std::sqrt(5.0);
155 
156  x1 = .942882415695480;
157  x2 = .641853342345781;
158  x3 = .236383199662150;
159  x(0) = a;
160  x(1) = m - x1 * h;
161  x(2) = m - alpha * h;
162  x(3) = m - x2 * h;
163  x(4) = m - beta * h;
164  x(5) = m - x3 * h;
165  x(6) = m;
166  x(7) = m + x3 * h;
167  x(8) = m + beta * h;
168  x(9) = m + x2 * h;
169  x(10) = m + alpha * h;
170  x(11) = m + x1 * h;
171  x(12) = b;
172 
173  y = apply_function<double>(f, x);
174 
175  fa = y(0);
176  fb = y(12);
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));
179 
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));
182 
183  s = sign(is);
184  if (s == 0.0)
185  s = 1;
186 
187  erri1 = std::abs(i1 - is);
188  erri2 = std::abs(i2 - is);
189 
190  R = 1;
191  if (erri2 != 0.0)
192  R = erri1 / erri2;
193 
194  if (R > 0 && R < 1)
195  tol2 = tol2 / R;
196 
197  is = s * std::abs(is) * tol2 / std::numeric_limits<double>::epsilon();
198  if (is == 0.0)
199  is = b - a;
200 
201  Q = quadlstep(f, a, b, fa, fb, is);
202 
203  return Q;
204 }
205 
206 
207 } // namespace itpp
SourceForge Logo

Generated on Fri Mar 21 2014 17:14:12 for IT++ by Doxygen 1.8.1.2