L'intégration numérique en C++

Autant que je sache, il n'existe pas de fonctions du type que vous cherchez dans la bibliothèque standard. Voici une mise en œuvre de l':

L'étendue de forme trapézoïdale à la règle.

Pour une fonction fixe f(x) à être intégré entre fixe
les limites de a et b, on peut doubler le nombre d'intervalles de la forme trapézoïdale de la règle, sans perdre le bénéfice du travail précédent. Les plus grossiers de la mise en œuvre de la trapézoïdale règle est en moyenne de la fonction au niveau de ses extrémités a et b. La première étape de raffinement est à ajouter à cette moyenne, la valeur de la fonction au point à mi-chemin. La deuxième étape de raffinement est d'ajouter de la valeur à la 1/4 et 3/4 points.

Un certain nombre de quadrature élémentaire algorithmes nécessitent l'ajout d'
les étapes successives de raffinement. Il est commode pour encapsuler cette fonctionnalité dans un Quadrature structure:

struct Quadrature
{  
   //Abstract base class for elementary quadrature algorithms.
   Int n; //Current level of refinement.

   virtual Doub next() = 0;
   //Returns the value of the integral at the nth stage of refinement. 
   //The function next() must be defined in the derived class.
};

Puis le Trapzd la structure est dérivée à partir de ce comme suit:

template<class T> 
struct Trapzd: Quadrature
{
Doub a, b, s; //Limits of integration and current value of integral.
T &func;
Trapzd() { };
//func is function or functor to be integrated between limits: a and b 
Trapzd(T &funcc, const Doub aa, const Doub bb)   
: func(funcc), a(aa), b(bb)
{
n = 0;
}
//Returns the nth stage of refinement of the extended trapezoidal rule. 
//On the first call (n = 1), the routine returns the crudest estimate  
//of integral of f x /dx in [a,b]. Subsequent calls set n=2,3,... and
//improve the accuracy by adding 2n - 2 additional interior points.
Doub next()
{
Doub x, tnm, sum, del;
Int it, j;
n++;
if (n == 1)
{
return (s = 0.5 * (b-a) * (func(a) + func(b)));
} 
else
{
for (it = 1, j = 1; j < n - 1; j++)
{
it <<= 1;
}
tnm = it;
//This is the spacing of the points to be added.          
del = (b - a) / tnm; 
x = a + 0.5 * del;
for (sum = 0.0,j = 0; j < it; j++, x += del)
{
sum += func(x);
}
//This replaces s by its refined value.  
s = 0.5 * (s + (b - a) * sum / tnm); 
return s;
}
}
};

Utilisation:

La Trapzd structure pourrait être utilisé de plusieurs façons. L'
la plus simple et la plus grossière est d'intégrer une fonction de l'étendue de forme trapézoïdale de la règle où vous savez à l'avance le nombre d'étapes que vous le souhaitez. Si vous
voulez 2^M + 1, vous pouvez accomplir cela en le fragment:

Ftor func;      //Functor func here has no parameters.
Trapzd<Ftor> s(func, a, b);
for(j = 1 ;j <= m + 1; j++) val = s.next();

avec la réponse retournée comme val. Ici Ftor est un foncteur contenant la fonction à intégrer.

Bonus:

Beaucoup mieux, bien sûr, est d'affiner la forme trapézoïdale de la règle jusqu'à un certain
degré de précision qui a été réalisé. Une fonction est la suivante:

template<class T>
Doub qtrap(T &func, const Doub a, const Doub b, const Doub eps = 1.0e-10)
{
//Returns the integral of the function or functor func from a to b. 
//The constants EPS can be set to the desired fractional accuracy and    
//JMAX so that 2 to the power JMAX-1 is the maximum allowed number of   
//steps. Integration is performed by the trapezoidal rule.
const Int JMAX = 20;
Doub s, olds = 0.0; //Initial value of olds is arbitrary.
Trapzd<T> t(func, a, b);
for (Int j = 0; j < JMAX; j++) 
{
s = t.next();
if (j > 5) //Avoid spurious early convergence.
{
if (abs(s - olds) < eps * abs(olds) || (s == 0.0 && olds == 0.0)) 
{
return s;
}
}
olds = s;
}
throw("Too many steps in routine qtrap");
}

Typedefs.

typedef double Doub;    //64 - bit floating point
typedef int Int;        //32 - bit signed integer

OriginalL'auteur Ziezi

Vous pouvez consulter le coup de pouce de Quadrature et de la Différenciation de la Bibliothèque. En particulier, ils fournissent une version de la Règle du Trapèze:

https://www.boost.org/doc/libs/1_69_0/libs/math/doc/html/math_toolkit/trapezoidal.html

La Quadrature/Différenciation de la Bibliothèque est très bien écrit et compatible avec le C++ moderne, en ce qu'il est possible d'utiliser une expression lambda ou un objet de fonction pour la fonction à intégrer. J'ai été opérationnel très rapidement.

Voici un exemple d'approximation de pi, par l'approximation de l'intégrale des 4/(1 + x^2), x = 0 à x = 1, en mettant l'intégrande comme une expression lambda.

#include <boost/math/quadrature/trapezoidal.hpp>
#include <iostream>
using boost::math::quadrature::trapezoidal;
using std::cout;
using std::endl;
//Put inside a test function:
auto f = [](double x)
{
return 4.0 / (1.0 + x * x);
};
double appPi = trapezoidal(f, 0.0, 1.0);
double tol = 1e-6;
int max_refinements = 20;
double appPi2 = trapezoidal(f, 0.0, 1.0, tol, max_refinements);
cout << "Trapezoid Rule results for computing pi by integration" << endl;
cout << "a) with defaults, and b) with tol and max_refine set : " << endl;
cout << appPi << ", " << appPi2 << endl << endl;

J'ai donné deux exemples, l'un utilisant les paramètres par défaut pour la discrétisation de la plage de l'intégration et de convergence, et la seconde à l'aide de paramètres personnalisés.

Mes résultats (juste de faire une copie de la sortie de l'écran):

Trapezoid Rule results for computing pi by integration
a) with defaults, and b) with tol and max_refine set :
3.14159, 3.14159

OriginalL'auteur djhanson