Numerische integration in C++

Soweit ich weiß gibt es keine Funktionen der Art, die Sie suchen, in der standard-Bibliothek. Hier ist eine Implementierung des:

Der erweiterten trapezregel.

Für eine Feste Funktion f(x) integriert werden zwischen fixen
Grenzen a und b eine Verdoppelung der Anzahl der Intervalle in der erweiterten trapezregel, ohne dabei den nutzen der bisherigen Arbeit. Das gröbste Umsetzung der Trapez-Regel ist der Mittelwert der Funktion an den Endpunkten a und b. Die erste Stufe der Verfeinerung ist, um zu diesem Mittelwert wird der Wert der Funktion im Punkt auf halbem Weg. Die zweite Stufe der Verfeinerung ist, fügen Sie die Werte in die 1/4 und 3/4 Punkte.

Einer Reihe von elementaren Quadratur-algorithmen beinhalten das hinzufügen
aufeinanderfolgenden Stufen der Verfeinerung. Es ist bequem zu Kapseln diese Funktion bei Quadrature Struktur:

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.
};

Dann die Trapzd Struktur leitet sich von diesem wie folgt:

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;
        }
    }
};

Verwendung:

Den Trapzd Struktur könnte auf mehrere Arten verwendet werden. Die
einfachste und rohsten zu integrieren, ist eine Funktion, mit der erweiterten trapezregel, wo Sie im Voraus wissen, die Anzahl der Schritte, die Sie wollen. Wenn Sie
wollen 2^M + 1 Sie können dies erreichen, indem das 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();

mit der Antwort zurück, wie val. Hier Ftor ist ein Funktor mit der Funktion integriert werden.

Bonus:

Viel besser, natürlich, ist die Verfeinerung der trapezregel bis auf einige angegeben
Grad der Genauigkeit erreicht wurde. Eine Funktion dafür ist:

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

InformationsquelleAutor Ziezi

Möchten Sie vielleicht zu prüfen, die Boost-Quadratur und Differenzierung Bibliothek. Insbesondere bieten Sie eine version der Trapez-Regel:

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

Die Quadratur/Differenzierung Bibliothek ist sehr gut geschrieben und kompatibel mit modernen C++, das kann man nur mit einem lambda-Ausdruck oder ein function-Objekt für den integranden um. Ich war und läuft sehr schnell.

Hier ist ein Beispiel der approximation pi, die Annäherung an das integral von 4/(1 + x^2), von x = 0 bis x = 1, setzen Sie den integranden um, da ein lambda-Ausdruck.

#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;

Ich zwei Beispiele, eines mit den Standardeinstellungen für die Diskretisierung der Bereich der integration und Konvergenz, und die zweite mit benutzerdefinierten Einstellungen.

Meine Ergebnisse (nur eine Kopie der Ausgabe auf dem Bildschirm) sind:

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

InformationsquelleAutor djhanson