Theoretica
A C++ numerical and automatic mathematical library
Loading...
Searching...
No Matches
integral.h
Go to the documentation of this file.
1
5
6#ifndef THEORETICA_INTEGRATION_H
7#define THEORETICA_INTEGRATION_H
8
9#include "../core/constants.h"
10#include "../core/function.h"
11#include "../polynomial/polynomial.h"
12#include "../polynomial/orthogonal.h"
13#include "./gauss.h"
14
15
16namespace theoretica {
17
18
23 template<typename Field = real>
25
27 P.coeff.resize(p.size() + 1);
28 P[0] = Field(0.0);
29
30 for (unsigned int i = 0; i < p.size(); ++i)
31 P[i + 1] = p[i] / Field(i + 1);
32
33 return P;
34 }
35
36
45 inline real integral(const polynomial<real>& p, real a, real b) {
46
47 real P_a = 0.0;
48 real P_b = 0.0;
49
50 for (unsigned int i = 0; i < p.size(); ++i) {
51
52 const unsigned int pos = p.size() - i - 1;
53 P_a = a * (p[pos] / (pos + 1) + P_a);
54 P_b = b * (p[pos] / (pos + 1) + P_b);
55 }
56
57 return P_b - P_a;
58 }
59
60
69 template<typename RealFunction>
71 unsigned int steps = CALCULUS_INTEGRAL_STEPS) {
72
73 if(steps == 0) {
74 TH_MATH_ERROR("integral_midpoint", steps, DIV_BY_ZERO);
75 return nan();
76 }
77
78 const real dx = (b - a) / steps;
79 real res = 0;
80
81 for (unsigned int i = 0; i < steps; ++i)
82 res += f(a + (i + 0.5) * dx);
83
84 return res * dx;
85 }
86
87
96 template<typename RealFunction>
98 unsigned int steps = CALCULUS_INTEGRAL_STEPS) {
99
100 if(steps == 0) {
101 TH_MATH_ERROR("integral_trapezoid", steps, DIV_BY_ZERO);
102 return nan();
103 }
104
105 const real dx = (b - a) / steps;
106 real res = 0;
107
108 res += 0.5 * (f(a) + f(b));
109
110 for (unsigned int i = 1; i < steps; ++i)
111 res += f(a + i * dx);
112
113 return res * dx;
114 }
115
116
125 template<typename RealFunction>
127 unsigned int steps = CALCULUS_INTEGRAL_STEPS) {
128
129 if(steps == 0) {
130 TH_MATH_ERROR("integral_simpson", steps, DIV_BY_ZERO);
131 return nan();
132 }
133
134 const real dx = (b - a) / (real) steps;
135 real res = 0;
136
137 // Sum terms by order of magnitude supposing that
138 // f stays at the same order inside the interval,
139 // to alleviate truncation errors
140
141 res += f(a) + f(b);
142
143 for (unsigned int i = 2; i < steps; i += 2)
144 res += 2 * f(a + i * dx);
145
146 for (unsigned int i = 1; i < steps; i += 2)
147 res += 4 * f(a + i * dx);
148
149 return res * dx / 3.0;
150 }
151
152
161 template<typename RealFunction>
163 RealFunction f,
164 real a, real b,
165 unsigned int iter = 8) {
166
167 real T[iter][iter];
168
169 T[0][0] = (f(a) + f(b)) * (b - a) / 2.0;
170
171 for (unsigned int j = 1; j < iter; ++j) {
172
173 // Composite Trapezoidal Rule
174 T[j][0] = integral_trapezoid(f, a, b, 1 << j);
175
176 // Richardson extrapolation
177 for (unsigned int k = 1; k <= j; ++k) {
178 const uint64_t coeff = uint64_t(1) << (2 * k);
179 T[j][k] = (coeff * T[j][k - 1] - T[j - 1][k - 1]) / (coeff - 1);
180 }
181 }
182
183 // Return the best approximation
184 return T[iter - 1][iter - 1];
185 }
186
187
196 template<typename RealFunction>
198 RealFunction f,
199 real a, real b,
200 real tolerance = CALCULUS_INTEGRAL_TOL) {
201
202 const unsigned int MAX_ROMBERG_ITER = 16;
204
205 T[0][0] = (f(a) + f(b)) * (b - a) / 2.0;
206
207 for (unsigned int j = 1; j < 16; ++j) {
208
209 // Composite Trapezoidal Rule
210 T[j][0] = integral_trapezoid(f, a, b, 1 << j);
211
212 // Richardson extrapolation
213 for (unsigned int k = 1; k <= j; ++k) {
214 const uint64_t coeff = uint64_t(1) << (2 * k);
215 T[j][k] = (coeff * T[j][k - 1] - T[j - 1][k - 1]) / (coeff - 1);
216 }
217
218 // Stop the algorithm when the desired precision has been reached
219 if(abs(T[j][j] - T[j - 1][j - 1]) < tolerance)
220 return T[j][j];
221 }
222
223 // Return the best approximation
224 return T[MAX_ROMBERG_ITER - 1][MAX_ROMBERG_ITER - 1];
225 }
226
227
233 template<typename RealFunction>
235 RealFunction f, const std::vector<real>& x, const std::vector<real>& w) {
236
237 if(x.size() != w.size()) {
238 TH_MATH_ERROR("integral_gauss", x.size(), INVALID_ARGUMENT);
239 return nan();
240 }
241
242 real res = 0;
243
244 for (int i = 0; i < x.size(); ++i)
245 res += w[i] * f(x[i]);
246
247 return res;
248 }
249
250
257 template<typename RealFunction>
259 RealFunction f, real* x, real* w, unsigned int n) {
260
261 real res = 0;
262
263 for (unsigned int i = 0; i < n; ++i)
264 res += w[i] * f(x[i]);
265
266 return res;
267 }
268
269
277 template<typename RealFunction>
279 RealFunction f, real* x, real* w, unsigned int n,
281
282 real res = 0;
283
284 for (unsigned int i = 0; i < n; ++i)
285 res += w[i] * f(x[i]) * Winv(x[i]);
286
287 return res;
288 }
289
290
300 template<typename RealFunction>
302 RealFunction f, real a, real b, real* x, real* w, unsigned int n) {
303
304 const real mean = (b + a) / 2.0;
305 const real halfdiff = (b - a) / 2.0;
306
307 real res = 0;
308
309 for (int i = n - 1; i >= 0; --i)
310 res += w[i] * f(halfdiff * x[i] + mean);
311
312 return res * halfdiff;
313 }
314
315
325 template<typename RealFunction>
327 RealFunction f, real a, real b,
328 const std::vector<real>& x, const std::vector<real>& w) {
329
330 const real mean = (b + a) / 2.0;
331 const real halfdiff = (b - a) / 2.0;
332
333 real res = 0;
334
335 for (int i = x.size() - 1; i >= 0; --i)
336 res += w[i] * f(halfdiff * x[i] + mean);
337
338 return res * halfdiff;
339 }
340
341
350 template<typename RealFunction>
352 RealFunction f, real a, real b, const std::vector<real>& x) {
353
354 return integral_legendre(f, a, b, x, legendre_weights(x));
355 }
356
357
368 template<typename RealFunction>
369 inline real integral_legendre(RealFunction f, real a, real b, unsigned int n = 16) {
370
371 switch(n) {
372 case 2: return integral_legendre(f, a, b,
373 tables::legendre_roots_2, tables::legendre_weights_2, 2); break;
374 case 4: return integral_legendre(f, a, b,
375 tables::legendre_roots_4, tables::legendre_weights_4, 4); break;
376 case 8: return integral_legendre(f, a, b,
377 tables::legendre_roots_8, tables::legendre_weights_8, 8); break;
378 case 16: return integral_legendre(f, a, b,
379 tables::legendre_roots_16, tables::legendre_weights_16, 16); break;
380 // case 32: return integral_legendre(f, a, b,
381 // tables::legendre_roots_32, tables::legendre_weights_32, 32); break;
382 default: return integral_legendre(f, a, b, legendre_roots(n)); break;
383 }
384 }
385
386
393 template<typename RealFunction>
394 inline real integral_laguerre(RealFunction f, const std::vector<real>& x) {
395
396 return integral_gauss(f, x, laguerre_weights(x));
397 }
398
399
408 template<typename RealFunction>
410 RealFunction f, real a, real b, const std::vector<real>& x) {
411
412 const std::vector<real> weights = laguerre_weights(x);
413
414 const real exp_a = exp(-a);
415 const real exp_b = exp(-b);
416
417 real res = 0;
418
419 for (int i = x.size() - 1; i >= 0; --i)
420 res += weights[i] * (exp_a * f(x[i] + a) - exp_b * f(x[i] + b));
421
422 return res;
423 }
424
425
434 template<typename RealFunction>
435 inline real integral_laguerre(RealFunction f, unsigned int n = 16) {
436
437 switch(n) {
438 case 2: return integral_gauss(f,
439 tables::laguerre_roots_2, tables::laguerre_weights_2, 2); break;
440 case 4: return integral_gauss(f,
441 tables::laguerre_roots_4, tables::laguerre_weights_4, 4); break;
442 case 8: return integral_gauss(f,
443 tables::laguerre_roots_8, tables::laguerre_weights_8, 8); break;
444 case 16: return integral_gauss(f,
445 tables::laguerre_roots_16, tables::laguerre_weights_16, 16); break;
446 // case 32: return integral_gauss(f,
447 // tables::laguerre_roots_32, tables::laguerre_weights_32, 32); break;
448 default: {
449 TH_MATH_ERROR("integral_laguerre", n, INVALID_ARGUMENT);
450 return nan(); break;
451 }
452 }
453 }
454
455
462 template<typename RealFunction>
463 inline real integral_hermite(RealFunction f, const std::vector<real>& x) {
464
465 return integral_gauss(f, x, hermite_weights(x));
466 }
467
468
477 template<typename RealFunction>
478 inline real integral_hermite(RealFunction f, unsigned int n = 16) {
479
480 switch(n) {
481 case 2: return integral_gauss(f,
482 tables::hermite_roots_2, tables::hermite_weights_2, 2); break;
483 case 4: return integral_gauss(f,
484 tables::hermite_roots_4, tables::hermite_weights_4, 4); break;
485 case 8: return integral_gauss(f,
486 tables::hermite_roots_8, tables::hermite_weights_8, 8); break;
487 case 16: return integral_gauss(f,
488 tables::hermite_roots_16, tables::hermite_weights_16, 16); break;
489 default: {
490 TH_MATH_ERROR("integral_hermite", n, INVALID_ARGUMENT);
491 return nan(); break;
492 }
493 }
494 }
495
496
502 real_function f, real a, real step_sz = 1,
503 real tol = CALCULUS_INTEGRAL_TOL, unsigned int max_iter = 100) {
504
505 // Current lower extreme of the interval
506 real x_n = a + step_sz;
507
508 // Total integral sum
510
511 // Variation between steps
512 real delta = inf();
513
514 // Number of steps performed
515 unsigned int i = 0;
516
517 while(abs(delta) > tol && i <= max_iter) {
518
520 sum += delta;
521 x_n += step_sz;
522 i++;
523 }
524
525 if(i >= max_iter) {
526 TH_MATH_ERROR("integral_inf_riemann", i, NO_ALGO_CONVERGENCE);
527 return nan();
528 }
529
530 return sum;
531 }
532
533
541 template<typename RealFunction>
543 return integral_romberg_tol(f, a, b);
544 }
545
546
556 template<typename RealFunction>
558 return integral_romberg_tol(f, a, b, tol);
559 }
560
561}
562
563
564#endif
#define TH_MATH_ERROR(F_NAME, VALUE, EXCEPTION)
TH_MATH_ERROR is a macro which throws exceptions or modifies errno (depending on which compiling opti...
Definition error.h:219
Roots and weights for Gaussian quadrature.
Main namespace of the library which contains all functions and objects.
Definition algebra.h:27
std::vector< real > hermite_weights(const std::vector< real > &roots)
Hermite weights for Gauss-Hermite quadrature of n-th order.
Definition orthogonal.h:316
double real
A real number, defined as a floating point type.
Definition constants.h:198
real inf()
Return positive infinity in floating point representation.
Definition error.h:76
dual2 abs(dual2 x)
Compute the absolute value of a second order dual number.
Definition dual2_functions.h:198
std::remove_reference_t< decltype(std::declval< Structure >()[0])> vector_element_t
Extract the type of a vector (or any indexable container) from its operator[].
Definition core_traits.h:134
dual2 exp(dual2 x)
Compute the exponential of a second order dual number.
Definition dual2_functions.h:138
std::vector< real > legendre_weights(const std::vector< real > &roots)
Legendre weights for Gauss-Legendre quadrature of n-th order.
Definition orthogonal.h:274
std::vector< real > laguerre_weights(const std::vector< real > &roots)
Laguerre weights for Gauss-Laguerre quadrature of n-th order.
Definition orthogonal.h:295
auto sum(const Vector &X)
Compute the sum of a vector of real values using pairwise summation to reduce round-off error.
Definition dataset.h:219
real integral_inf_riemann(real_function f, real a, real step_sz=1, real tol=CALCULUS_INTEGRAL_TOL, unsigned int max_iter=100)
Integrate a function from a point up to infinity by integrating it by steps, stopping execution when ...
Definition integral.h:501
real integral_romberg(RealFunction f, real a, real b, unsigned int iter=8)
Approximate the definite integral of an arbitrary function using Romberg's method accurate to the giv...
Definition integral.h:162
constexpr int CALCULUS_INTEGRAL_STEPS
Default number of steps for integral approximation.
Definition constants.h:282
real integral_gauss(RealFunction f, const std::vector< real > &x, const std::vector< real > &w)
Use Gaussian quadrature using the given points and weights.
Definition integral.h:234
real integral_laguerre(RealFunction f, const std::vector< real > &x)
Use Gauss-Laguerre quadrature of arbitrary degree to approximate an integral over [0,...
Definition integral.h:394
real integral_trapezoid(RealFunction f, real a, real b, unsigned int steps=CALCULUS_INTEGRAL_STEPS)
Approximate the definite integral of an arbitrary function using the trapezoid method.
Definition integral.h:97
real integral_romberg_tol(RealFunction f, real a, real b, real tolerance=CALCULUS_INTEGRAL_TOL)
Approximate the definite integral of an arbitrary function using Romberg's method to the given tolera...
Definition integral.h:197
polynomial< Field > integral(const polynomial< Field > &p)
Compute the indefinite integral of a polynomial.
Definition integral.h:24
real integral_hermite(RealFunction f, const std::vector< real > &x)
Use Gauss-Hermite quadrature of arbitrary degree to approximate an integral over (-inf,...
Definition integral.h:463
real integral_simpson(RealFunction f, real a, real b, unsigned int steps=CALCULUS_INTEGRAL_STEPS)
Approximate the definite integral of an arbitrary function using Simpson's method.
Definition integral.h:126
real integral_legendre(RealFunction f, real a, real b, real *x, real *w, unsigned int n)
Use Gauss-Legendre quadrature of arbitrary degree to approximate a definite integral providing the ro...
Definition integral.h:301
std::function< real(real)> real_function
Function pointer to a real function of real variable.
Definition function.h:20
real integral_midpoint(RealFunction f, real a, real b, unsigned int steps=CALCULUS_INTEGRAL_STEPS)
Approximate the definite integral of an arbitrary function using the midpoint method.
Definition integral.h:70
real nan()
Return a quiet NaN number in floating point representation.
Definition error.h:54
std::vector< real > legendre_roots(unsigned int n)
Roots of the n-th Legendre polynomial.
Definition orthogonal.h:249