6#ifndef THEORETICA_INTEGRATION_H
7#define THEORETICA_INTEGRATION_H
9#include "../core/constants.h"
10#include "../core/function.h"
11#include "../polynomial/polynomial.h"
12#include "../polynomial/orthogonal.h"
23 template<
typename Field = real>
27 P.coeff.resize(
p.size() + 1);
30 for (
unsigned int i = 0; i <
p.size(); ++i)
50 for (
unsigned int i = 0; i <
p.size(); ++i) {
52 const unsigned int pos =
p.size() - i - 1;
69 template<
typename RealFunction>
81 for (
unsigned int i = 0; i <
steps; ++i)
82 res += f(a + (i + 0.5) *
dx);
96 template<
typename RealFunction>
108 res += 0.5 * (f(a) + f(b));
110 for (
unsigned int i = 1; i <
steps; ++i)
111 res += f(a + i *
dx);
125 template<
typename RealFunction>
143 for (
unsigned int i = 2; i <
steps; i += 2)
144 res += 2 * f(a + i *
dx);
146 for (
unsigned int i = 1; i <
steps; i += 2)
147 res += 4 * f(a + i *
dx);
149 return res *
dx / 3.0;
161 template<
typename RealFunction>
165 unsigned int iter = 8) {
169 T[0][0] = (f(a) + f(b)) * (b - a) / 2.0;
171 for (
unsigned int j = 1;
j <
iter; ++
j) {
177 for (
unsigned int k = 1;
k <=
j; ++
k) {
179 T[
j][
k] = (coeff *
T[
j][
k - 1] -
T[
j - 1][
k - 1]) / (coeff - 1);
196 template<
typename RealFunction>
205 T[0][0] = (f(a) + f(b)) * (b - a) / 2.0;
207 for (
unsigned int j = 1;
j < 16; ++
j) {
213 for (
unsigned int k = 1;
k <=
j; ++
k) {
215 T[
j][
k] = (coeff *
T[
j][
k - 1] -
T[
j - 1][
k - 1]) / (coeff - 1);
233 template<
typename RealFunction>
235 RealFunction f,
const std::vector<real>& x,
const std::vector<real>&
w) {
237 if(x.size() !=
w.size()) {
244 for (
int i = 0; i < x.size(); ++i)
245 res +=
w[i] * f(x[i]);
257 template<
typename RealFunction>
263 for (
unsigned int i = 0; i <
n; ++i)
264 res +=
w[i] * f(x[i]);
277 template<
typename RealFunction>
284 for (
unsigned int i = 0; i <
n; ++i)
300 template<
typename RealFunction>
304 const real mean = (b + a) / 2.0;
309 for (
int i =
n - 1; i >= 0; --i)
325 template<
typename RealFunction>
328 const std::vector<real>& x,
const std::vector<real>&
w) {
330 const real mean = (b + a) / 2.0;
335 for (
int i = x.size() - 1; i >= 0; --i)
350 template<
typename RealFunction>
368 template<
typename RealFunction>
373 tables::legendre_roots_2, tables::legendre_weights_2, 2);
break;
375 tables::legendre_roots_4, tables::legendre_weights_4, 4);
break;
377 tables::legendre_roots_8, tables::legendre_weights_8, 8);
break;
379 tables::legendre_roots_16, tables::legendre_weights_16, 16);
break;
393 template<
typename RealFunction>
408 template<
typename RealFunction>
419 for (
int i = x.size() - 1; i >= 0; --i)
434 template<
typename RealFunction>
439 tables::laguerre_roots_2, tables::laguerre_weights_2, 2);
break;
441 tables::laguerre_roots_4, tables::laguerre_weights_4, 4);
break;
443 tables::laguerre_roots_8, tables::laguerre_weights_8, 8);
break;
445 tables::laguerre_roots_16, tables::laguerre_weights_16, 16);
break;
462 template<
typename RealFunction>
477 template<
typename RealFunction>
482 tables::hermite_roots_2, tables::hermite_weights_2, 2);
break;
484 tables::hermite_roots_4, tables::hermite_weights_4, 4);
break;
486 tables::hermite_roots_8, tables::hermite_weights_8, 8);
break;
488 tables::hermite_roots_16, tables::hermite_weights_16, 16);
break;
526 TH_MATH_ERROR(
"integral_inf_riemann", i, NO_ALGO_CONVERGENCE);
541 template<
typename RealFunction>
556 template<
typename RealFunction>
#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