Numerische Integration mit boost?



  • Hallo,

    ich suche eine Möglichkeit mit boost eine eindimensionale Funktionen numerisch zu integrieren. Das einzige, was ich bisher dazu gefunden habe ist boost::numeric::odeint. Obwohl es sehr viele Beispiele dazu gibt, sieht es für mich so aus, als könnte man damit nur ODEs lösen. Wie ich jedoch ganz normale Funktionen damit integrieren kann, erschließt sich mir jedoch nicht.
    Kann mir jemand erklären, was es genauer damit auf sich hat? 🙂



  • Sodele schrieb:

    Hallo,

    ich suche eine Möglichkeit mit boost eine eindimensionale Funktionen numerisch zu integrieren. Das einzige, was ich bisher dazu gefunden habe ist boost::numeric::odeint. Obwohl es sehr viele Beispiele dazu gibt, sieht es für mich so aus, als könnte man damit nur ODEs lösen. Wie ich jedoch ganz normale Funktionen damit integrieren kann, erschließt sich mir jedoch nicht.
    Kann mir jemand erklären, was es genauer damit auf sich hat? 🙂

    Du verwendest einfach nicht den Zustandsvektor sondern den "Zeit"-Wert zur Berechnung. "Normale" Funktionen sind unabhängig vom vorherigen Zustand, es zählt nur der aktuelle "Zeit"-Schritt. Zum Beispiel so für die Funktion

    f(x)=x2f(x) = x^2

    Wenn du die Funktion numerisch von (beispielsweise) 0 bis 10 integrieren willst, geht das so:

    #include <iostream>
    #include <array>
    #include <boost/numeric/odeint.hpp>
    
    using namespace boost::numeric::odeint;
    using state_type = std::array < double, 1 > ;
    
    void func(const state_type &x, state_type &dxdt, double t)
    {
    	dxdt[0] = t*t;
    }
    
    void observer(const state_type &x, const double t)
    {
    	std::cout << t << '\t' << x[0] << '\n';
    }
    
    int main()
    {
    	state_type x = { 0 };
    	integrate(func, x, 0.0, 10.0, 0.1, observer);
    }
    

    Das Ergebnis ist 10003333,333\frac{1000}{3}\approx 333,333, wie das Program auch ausgibt und dieses ist nach dem Aufruf von integrate auch in der Variable x gespeichert. Wenn du nur das Ergebnis brauchst kannst du den Observer natürlich auch weglassen.



  • Danke, jetzt habe auch ich es verstanden.
    Eine Frage noch, wie soll ich den Parameter dt am besten Wählen? In der boost-Dokumentation steht folgendes:

    Integration start at t0 and x0 and ends at some t' = t0 + n dt with n such that t1 - dt < t' <= t1.

    Jedoch kann ich daraus nicht schließen, wie groß man dt für ein möglichst optimales Ergebniss wählen muss.



  • Sodele schrieb:

    Danke, jetzt habe auch ich es verstanden.
    Eine Frage noch, wie soll ich den Parameter dt am besten Wählen?

    Viele der stepper aus odeint haben eine adaptive Schrittweite, d.h. du musst dich (fast) nicht weiter drum kümmern. In dem geposteten Beispiel wird auch ein solcher verwendet (Dormand-Prince, auch der default Stepper in Matlab).

    Die Schrittweite sollte halt ungefähr passen (heißt: auf jeden Fall kleiner als das zu integrierende Interval), mit 0.1 fährt man meistens sehr gut. Falls die nicht passt, wird sie dann durch den stepper eben angepasst, wie man auch in an dem Output des observers sehen kann.

    Eine allgemein gültige Lösung die immer optimal ist gibt es für dieses Problem nicht. Ist aber auch nicht so schlimm, weil sich eigentlich nur die Laufzeit ein bisschen erhöht wenn deine Anfangsschrittweite sehr weit weg von der optimalen Schrittweite ist.


Log in to reply