Numerische Bibliothek in C++

headmyshoulder

hustbaer schrieb:

Hab mir interessehalber den Code etwas angesehen, dabei ist mir was aufgefallen.

// http://yanl.svn.sourceforge.net/viewvc/yanl/trunk/yanl/odeint/solver_runge_kutta4.hpp?revision=19&view=markup :
         std::vector<T> dxdt , dxt , dxm , xt ;

         void resize( size_t new_size ) {
             dxdt.resize( new_size , T(0.0) );
             dxt.resize( new_size , T(0.0)  ); // angenommen hier
             dxm.resize( new_size , T(0.0)  ); // oder hier
             xt.resize( new_size , T(0.0)  ); // oder hier fliegt eine exception
         }

         void next_step( CallMethod deriv, T t , T dt , T* state , size_t n) {
             if( n != dxdt.size() ) resize( n );
             // dann gibt's hier ein problem.
             //    -> lieber xt.size() abfragen
             ...
        }

Was genau meinst Du? Falls da ein exception fliegt crasht alles, solange man von aussen nichts fängt, oder irre ich mich?

hustbaer schrieb:

BTW: wieso ist der code so unübersichtlich formatiert?

Über Codeformatierung haben wir uns noch keine Gedanken gemacht, aber so schlecht sieht das nicht aus, oder?

headmyshoulder

Walli schrieb:

Hast Du die mal ausprobiert? Besonders die Sparse-Implementierung interessiert mich. Ich mache manchmal serielle Testprogrämmchen und benutze wegen der Performance immer einen selbstgefrickelten PETSc-Wrapper weil uBLAS grottenlangsam ist und MTL4 noch nicht wirklich stable.

@OP: Wenn Du wirklich viele Leute (inkl. mich) glücklich machen willst, dann schreibe einen gescheiten C++-Wrapper für PETSc oder eine saubere C++-Matrix-Lib, die im Prinzip mit geringem Aufwand alles mögliche interfacen kann . Ich wollte das immer schon mal ernsthaft angehen, aber irgendwie finde ich die Zeit nicht.

http://www.mcs.anl.gov/petsc/petsc-as/

Die Hardcore Lineare Algebra wollen wir nicht unbedingt bearbeiten. Das ist echt sehr advanced und man braucht ne Menge Erfahrung, was Compiler können und muss die Mathematik 100%ig verstehen. Obwohl es natürlich sehr interessant ist.

rüdiger

@Walli
ich hab Eigen bisher nur für Kleinigkeiten benutzt. Daher kann ich nicht wirklich sagen, ob es alles hält was es verspricht. Aber bisher finde ich es ganz gut.

@headmyshoulder
Willst du das Projekt eher zum lernen machen oder für eine wirkliche Aufgabe?

das solltet ihr euch auf jeden Fall anschauen
http://www.nr.com/

Die ältere 2. Edition gibt es auch kostenlos: http://www.nrbook.com/a/bookcpdf.php

hustbaer

headmyshoulder schrieb:

hustbaer schrieb:

Hab mir interessehalber den Code etwas angesehen, dabei ist mir was aufgefallen.

// http://yanl.svn.sourceforge.net/viewvc/yanl/trunk/yanl/odeint/solver_runge_kutta4.hpp?revision=19&view=markup :
         std::vector<T> dxdt , dxt , dxm , xt ;
 
         void resize( size_t new_size ) {
             dxdt.resize( new_size , T(0.0) );
             dxt.resize( new_size , T(0.0)  ); // angenommen hier
             dxm.resize( new_size , T(0.0)  ); // oder hier
             xt.resize( new_size , T(0.0)  ); // oder hier fliegt eine exception
         }

         void next_step( CallMethod deriv, T t , T dt , T* state , size_t n) {
             if( n != dxdt.size() ) resize( n );
             // dann gibt's hier ein problem.
             //    -> lieber xt.size() abfragen
             ...
        }

Was genau meinst Du? Falls da ein exception fliegt crasht alles, solange man von aussen nichts fängt, oder irre ich mich?

Klar crasht alles wenn eine Exception fliegt, die keiner fängt.
Aber es könnte sie ja auch jemand fangen.

Ich meine folgende Situation:

next_step() wird aufgerufen, und z.B. dxt.resize() schmeisst eine Exception.
die Exception wird gefangen
next_step() wird nochmal aufgerufen - mit dem gleichen Wert für "n" wie beim ersten Aufruf
next_step() vergleicht dann die Grösse von dxdt mit n, und sieht dass sie gleich sind. Und ruft daher resize() nicht nochmals auf.

Wenn "n" nun grösser ist, als die aktuelle Grösse von dxt, dmx oder xt, dann crasht der next_step() wahrscheinlich mit einer Access-Violation statt einer hübschen Exception. Oder noch schlimmer: liefert falsche Ergebnisse.

BTW: die "Lösung" die ich vorgeschlagen habe, ist auch nicht wirklich ausreichend.
Noch besser so in der Art:

template< class T, class CallMethod=dynamical_system<T>& >
     class ode_solver_runge_kutta4 : public ode_solver_base<T, CallMethod>
     {
         std::vector<T> m_dxdt;
         std::vector<T> m_dxt;
         std::vector<T> m_dxm;
         std::vector<T> m_xt;
         size_t m_buffer_size;

         void resize(size_t new_size)
         {
             if (new_size != m_buffer_size)
             {
                 m_buffer_size = 0; // if we're interrupted, force a re-start next time we're called
                 m_dxdt.resize(new_size, T(0.0));
                 m_dxt.resize(new_size, T(0.0));
                 m_dxm.resize(new_size, T(0.0));
                 m_xt.resize(new_size, T(0.0));
                 m_buffer_size = new_size;
             }
         }

     public:

         ode_solver_runge_kutta4() :
             m_buffer_size(0)
         {
         }

         void next_step(CallMethod deriv, T t , T dt , T* state , size_t n)
         {
             resize(n);

             T dh = dt * T(0.5)
             T d6 = dt / T(6.0);
             T th = dt * T(0.5) + t;

             deriv(state , m_dxdt.data() , t);
             for (size_t i = 0 ; i < n ; i++)
                 m_xt[i] = state[i] + dh * m_dxdt[i];

             deriv(m_xt.data() ,m_dxt.data(), th);
             for (size_t i = 0 ; i < n ; i++)
                 m_xt[i] = state[i] + dh * m_dxt[i];

             deriv(m_xt.data(), m_dxm.data(), th);
             for (size_t i = 0 ; i < n ; i++)
             {
                 m_xt[i] = state[i] + dt * m_dxm[i];
                 m_dxm[i] += m_dxt[i];
             }

             deriv(m_xt.data(), m_dxt.data(), t + dt);
             for (size_t i = 0; i < n ; i++)
                 state[i] += d6 * (m_dxdt[i] + m_dxt[i] + T(2.0) * m_dxm[i]);
         }
     };

headmyshoulder

hustbaer schrieb:

Ich meine folgende Situation:

next_step() wird aufgerufen, und z.B. dxt.resize() schmeisst eine Exception.

die Exception wird gefangen

next_step() wird nochmal aufgerufen - mit dem gleichen Wert für "n" wie beim ersten Aufruf
next_step() vergleicht dann die Grösse von dxdt mit n, und sieht dass sie gleich sind. Und ruft daher resize() nicht nochmals auf.

Ja, du hast Recht, das muß man beachten und ändern. Ich glaub dein Vorschlag ist auch eleganter als nach der Größe von xt zu fragen. Das wird geändert.

Ich glaub ich weiss jetzt auch was Du mit unschöner Codeformatierung meinst, obwohl ich members mit m_ zu bezeichnen nicht so mag.

headmyshoulder

rüdiger schrieb:

@headmyshoulder
Willst du das Projekt eher zum lernen machen oder für eine wirkliche Aufgabe?

Ich sehe dieses Projekt als richtige Aufgabe, wo was vernünftiges und benutzbares rauskommen soll. Die Idee ist ja nicht aus Spass an der Freude gekommen, sondern weil wir solch eine Bibliothek oft benutzen werden und uns die existierenden Lösungen nicht sehr gut gefallen.

Aber natürlich wird man dabei auch eine Menge lernen, sowohl thematisch als auch programmiertechnisch.

rüdiger schrieb:

das solltet ihr euch auf jeden Fall anschauen
http://www.nr.com/

Die ältere 2. Edition gibt es auch kostenlos: http://www.nrbook.com/a/bookcpdf.php

Ja, die Numerical Reciepes sind das Standardbuch für solche Probleme. Ich kenn das fast auswendig^^, zumindest einige Kapitel.

headmyshoulder

Mal ein kleines Update zum Stand der Dinge. Wir haben uns jetzt erstmal auf einen Part der lib konzentriert - odeint zum Lösen von Differentialgleichungen. Und der befindet sich in der boost sandbox, dort wird jetzt auch entwickelt. Wir haben jetzt ca. 2 Wochen an 20 Zeilen Code diskutiert und 1000mal wieder neuangefangen, aber ich denke wir sind jetzt auf einem guten Weg.

Falls jemand Lust hat, daran mitzuarbeiten würden wir uns sehr freuen. Besonders suchen wir jemanden, der sich richtig gut mit templates und metaprogrammierung auskennt.

knivil

Lösen von Differentialgleichungen

Da habt ihr euch aber etwas einfaches herausgesucht, dass die meisten Physiker/Mathematiker sicher selbst irgendwann mal implementiert haben (neben FFT).

headmyshoulder

knivil schrieb:

Lösen von Differentialgleichungen

Da habt ihr euch aber etwas einfaches herausgesucht, dass die meisten Physiker/Mathematiker sicher selbst irgendwann mal implementiert haben (neben FFT).

Joa, aber RK78 hat nicht jeder und Schrittweitensteuerung schreibt sich auch nicht jeder. Eine FFT würd ich nicht selber schreiben, die ist dann wahrscheinlich 1/10tel langsamer als die standards aus der fftw.

knivil

RK78 ... Schrittweitensteuerung

Ich musste ein Numerikbeleg fuer das Vordiplom machen. Dabei waehlte ich die Realisierung von allgemeinen expliziten Verfahren. Man konnte sie ueber das Butcher spezifiziert werden. Somit war es auch ein leichtes, Verfahren beliebiger Ordnung zu nutzen. Einfach Koeffizienten eintragen/laden. Inklusive Schrittweitensteuerung, beliebig dimensional. Jemand anderes hatte diagonalimplizite und implizite Verfahren implementiert, auch ueber ein Butcher-Schema definierbar. Und das im Grundstudium. Zur Performance ... naja so gut, wie man mit std::vector und std::list halt kommt. Aber bei meinem Programm bestand die Moeglichkeit, dynamisch Bibliotheken nachzuladen, die vielleicht SSE nutzten oder aber ergaenzende Verfahren bereitstellten, die nicht ins Butcher-Schema passten. Gekoppelt an gnuplot fuer Graphik und durch Lua scriptable (gut, das ist nicht euer Ziel).

Was aber nicht jeder hat, sind symplektische Verfahren (mit Schrittweitensteuerung) oder aber die Kontrolle der Loesung durch den globalen Diskretisierungsfehler. Oder aber ein Verfahren/Mapping von einer Schrittweiten gesteuerten Loesung auf eine aequidistant eingeteilte Zeitachse.

Was ich damit sagen will: die Konkurrenz ist gross, falls es nicht nur ein Selbstlernprojekt sein soll. Wenn es doch nur ein Selbstlernprojekt ist, dann hat es aber nichts in boost zu suchen.

Wieviel Mathematiker oder Physiker habt ihr denn in eurem Team? Damit meine ich jetzt nicht Experimentalphysiker, sondern Leute, die es taeglich einsetzen (wuerden), z.B. nichtlineare Dynamik oder Numerik. Was ist eure Zielgruppe?

Zu Templates: wo sollen sie eingesetzt werden und wo liegen die Vorteile?

headmyshoulder

knivil schrieb:

Ich musste ein Numerikbeleg fuer das Vordiplom machen. Dabei waehlte ich die Realisierung von allgemeinen expliziten Verfahren. Man konnte sie ueber das Butcher spezifiziert werden. Somit war es auch ein leichtes, Verfahren beliebiger Ordnung zu nutzen. Einfach Koeffizienten eintragen/laden. Inklusive Schrittweitensteuerung, beliebig dimensional. Jemand anderes hatte diagonalimplizite und implizite Verfahren implementiert, auch ueber ein Butcher-Schema definierbar. Und das im Grundstudium. Zur Performance ... naja so gut, wie man mit std::vector und std::list halt kommt.

Das klingt doch gut. Wenn Du Lust hast kannst Du die expliziten Verfahren ja in dem odeint Framework verallgemeinern.

knivil schrieb:

Was aber nicht jeder hat, sind symplektische Verfahren (mit Schrittweitensteuerung) oder aber die Kontrolle der Loesung durch den globalen Diskretisierungsfehler. Oder aber ein Verfahren/Mapping von einer Schrittweiten gesteuerten Loesung auf eine aequidistant eingeteilte Zeitachse.

Soweit ich weiss gibt es keine allgemeinen symplektischen Verfahren, sondern die hängen vom speziellen Problem ab. Es gibt einige Algorithmen für Hamiltonsysteme die sich in H = F(P) + G(Q) aufspalten lassen. Dafür haben wir auch Algorithmen, die einfach zu implementieren sind.

headmyshoulder

Hab ich grad erst gesehen:

knivil schrieb:

RK78 ... Schrittweitensteuerung

Was ich damit sagen will: die Konkurrenz ist gross, falls es nicht nur ein Selbstlernprojekt sein soll. Wenn es doch nur ein Selbstlernprojekt ist, dann hat es aber nichts in boost zu suchen.

Wieviel Mathematiker oder Physiker habt ihr denn in eurem Team? Damit meine ich jetzt nicht Experimentalphysiker, sondern Leute, die es taeglich einsetzen (wuerden), z.B. nichtlineare Dynamik oder Numerik. Was ist eure Zielgruppe?

Also, wir sind 4 Phyisker, alle aus dem Bereich Nichlineare Dynamik. Ein Selbstlernprojekt ist das definitiv nicht. Was uns nervt, das alle libs ihre eigenen Standards verwenden (eigene array, vektor) Klassen verwenden und das soll mit diesem Projekt umgangen werden. Gerade Anfänger in C/C++ die Numerik machen, nehmen sich ein Lehrbuch, wo beispielsweise sehr schön vector<double> erklärt wird. Dann möchte er eine Numerische bibliothek benutzen und stellt fest, dass man vorher seine vector<double>'s in array_xyz konvertieren muss. Das soll hier anders sein. Die Zielgruppe sind also alle Leute die Numerik machen. Natürlich ist der Anspruch nicht, für alles ein Lösung zu schaffen, z.b. NavierStokes Solver o.Ä., aber die Standards sollen behandelt werden. Irgendwo in dem Thread hab ich das auch noch ausführlicher besprochen.

knivil schrieb:

Zu Templates: wo sollen sie eingesetzt werden und wo liegen die Vorteile?

Templates werden zur Abstraktion eingesetzt. Die Algorithm sollen auf Containern arbeiten, beispielsweise vector,list,array. Und templates sind sau schnell. Wie schneller als ein Klassendesign.

knivil

Templates werden zur Abstraktion eingesetzt.

Hmm, es gibt noch andere Wege ...

Ich habe mir mal euren Entwurf in http://www.c-plusplus.net/forum/viewtopic-var-t-is-252700-and-highlight-is-.html angesehen. Ich werde mal in diesem Thread antworten, da es eher ein Designproblem ist, als ein Templateproblem. Um Containerunabhaengig zu sein, ist der Weg ueber Templates der falseche Weg. Genau dafuer wurden Iteratoren geschaffen. Wahrscheinlich ist der Grund:

wir sind 4 Phyisker

Was mir sofort beim Eulerschema auffaelt: wie soll eine variable Schrittweite beruecksichtigt werden? Das ganze nochmal Neueschreiben, nur mit Schrittweitensteuerung ... ist eigentlich unnoetig. Feste Schrittweite ist nur ein Spezialfall von variabler Schrittweite. Aber fangen wir am anfang an ... Es folgt mehr oder weniger C++ Pseudocode.

Eine gewoehnliche Dgl sieht z.B. so aus: dy/dt = f(y,t). Wie kann man das modellieren? Z.B. als Funktion oder Objekt. Ich tendiere ja eher zum Objekt, da ueber einen Konstruktor die Dgl noch parametrisiert werden kann.

class Ode
{
  public:
    // returns dy/dt = f( y, t)
    virtual void f(const InputContainer& in, double time, OutputContainer& out) const = 0;
};

Wobei ich mich hier auf eine Containerklasse schon festgelegt haette. Wahlweise (wie oben beschrieben) kann fuer den InputContainer auch ein Begin- und End-Iterator uebergeben werden und fuer den Outputcontainer halt ein einzelner Iterator, der zum schreiben in einen Container benutzt wird. Hier sind auch die Anforederungen am geringsten, da ein einzelner Forwarditerator ausreicht (ich wuerde aber dennoch ein RandomAccessIterator spezifizieren). Diese Schnittstelle ist kann dann mit std::vector, list oder gar c-Arrays genutzt werden. Zu den Kosten einer virtuellen Methode bitte die die entsprechenden Stellen im Meyers nachschlagen. Ein Ode(system) erbt von dieser abstrakten Klasse und ueberschreibt f.

Als zweites benoetigen Benoetigen wir ein Verfahren um von y_n nach y_(n+1) zu gelangen. Ein Beispiel waere das Eulerschema. In Ermangelung eines besseren Namens, nenne ich es mal Step.

class Step
{
  public:
    // calculates y_(n+1) out of y_n, t and time step dt
    virtual void next( const Ode& ode, input, t, dt, output ) const = 0;
};

Aus input werden wieder zwei Iteratoren, output aehnlich wie oben, t und dt einfache doubles. Eine Implementation fuer Euler koennte dann wie folgt aussehen:

class EulerStep : public Step
{
  public:
    virtual void next( ... )
    {
      ...
    }
};

Vorteil: maximale Flexibilitaet. Aktuell arbeitet ein Kollege an einem Problem, bei dem (vereinfacht) im 3 dimensinalen Raum die exakte Loesung fuer x und y bekannt ist (also keine Naeherungsloesung noetig), aber fuer z halt nicht. So kann einfach von Step fuer diesen Spezialfall abgeleitet werden, x und y exakt berechnet werden und z mittels was anderem approximiert werden. Das kann dann einfach mit einem bestehenden Iterationsscheam verwendet werden.

class Iteration
{
  public:
    virtual void iterate( const Ode& ode, const Step& step, start, end ) const = 0
};

Wobei start und end einfach das Zeitintervall angeben.

Fuer Schrittweitensteuerung kann man sich noch ueberlegen, eine Schrittweitenfunktion/Objekt (class StepController) entweder bei der Konstruktion eines Iterationobjektes als Parameter mitzugeben, oder aber der Funktionsparameter in iterate. Ich wuerde die Construktorvariante waehlen. Gleiches gilt fuer minimalen/maximalen Fehler und aehnliches.

Will man nun ein konkretes Verfahren wie das Heun ohne Schrittweitensteuerung kapselt man das ganze wieder in einer Klasse. Wobei die Funktion zur Ermittlung der Schrittweite eben einfach eine Konstante auf die aktuelle Zeit addiert.

Mit diesem Baukasten kann man glaube alle Einschrittverfahren abdecken. Und ist gleichzeitig sehr flexibel. Mit den einzelnen Legosteinen setzt man sein Loesungsverfahren zusammen. Wird ein Problem nicht von den mitgelieferten Komponenten abgedeckt, so koennen leicht neue oder angepasste entworfen werden. Einziges Problem ist, dass der Iterator vorher auf einen Container spezialisiert werden muss. Alternativ kann jedoch ein nuer Iteratortyp mit einem Interface entworfen werden, deren abgeleitet Klassen dann auf einen bestimmten Container durch Vererbung spezialisiert werden.

That's the big picture. So haette ich es gemacht. Dieser bjektorientierte Ansatz hat natuerlich auch Kosten. Sie sind in den Aufruf von virtuellen Funktionen zu suchen. Normal ist das aber kein Problem bei Einfachvererbung.

headmyshoulder

knivil schrieb:

...

Oha, find ich ja gut, daß es gleich konkret wird. Danke für die Kritik, viele der Sachen die Du angespricht, haben wir auch schon diskutiert.

Geschwindigkeit: Schau Dir mal die Speedtests in https://yanl.svn.sourceforge.net/svnroot/yanl/trunk/examples/odeint an. Das Klassendesign mit virtuellen Funktionen ist Faktor 2 langsamer^^. Sry, falls der Code dort leicht unübersichtlich ist.

Zum Design: Keiner hindert Dich dadran trotzdem Klassen zu verwenden. Hauptsache der (,,) operator ist definiert. Zum Beispiel kannst Du das obige Beispiel auch mit einer Klasse schreiben:

class Lorenz
{
public:

    void operator()( state_type &x , state_type &dxdt , double t )
    {
	dxdt[0] = sigma * ( x[1] - x[0] );
	dxdt[1] = R * x[0] - x[1] - x[0] * x[2];
	dxdt[2] = x[0]*x[1] - b * x[2];
    }
};

und das sieht dann angewandt so aus:

ode_step_euler< vector<double> > euler;
Lorenz lorenzo;
euler.next_step( lorenzo , x , t , dt );

Das ist doppelt so flexibel als ein Klassenansatz, weil man Funktion und auch Klassen benutzen kann.

Iteratoren oder Container: Ich denke das es mindestens eine Variante mit Iteratoren geben wird. Das obige Beispiel in dem anderen Thread ist nicht die Endfassung. Die Diskussion ist noch nicht zu Ende, beide Varianten haben Vor- und Nachteile. Wir können diese Diskussion ja woanders hinverschieben. Der grosste Nachteil ist, daß man die Anzahl der Elemente zwischen zwei Iteratoren nicht notwendigerweise mit O(0) rauskriegt. Bei Listen dauerts beispielsweise O(n). Falls irgendwann der Stepper weiterbenutzt werden soll, muß dort dann halt ein template Parameter stehen, statt einer abstrakten Klasse. Der entsprechende Typ muss die Stepper Eigenschaften erfüllen, genauso wie das dynamische System den (,,) Operator haben muss.

Schrittweitensteuerung Es gibt noch gar keine Schrittweitensteuerung. Fakt ist aber, daß falls die Schrittweite gefunden wurde, der Schritt auch ausgeführt werden soll und dafür ist die Funktionalität da mit beispielsweise ode_step_euler . Die eigentliche Iteration soll dann auch so ähnlich wie Du vorgeschlagen hast aussehen, nur daß noch Zustandsbeobachter hingezugefügt werden, welche in jeder Iteration irgendwas mit dem Zustand machen, Ausgabe, Berechnungen, Mitteln...

knivil

Sieht interessant aus, vom Ansatz her. Ausprobiert habe ich den Code noch nicht, nur gelesen. Aber etwas zum Codestyle: Beginnt man nach einer geschweiften Klammer eine neue Zeile, so rueckt man ein. Das macht das Lesen angenehmer. Ihr koennt euch ja mal einige Syleguides durchlesen.

Auch solltet ihr euch um geeignete Matrix und Vektorklassen bemuehen, die solche Zeilen kapseln:

for( size_t i=0 ; i<dxdt.size() ; ++i ) state[i] += dt * dxdt[i];

und in etwas Vertrauteres verwandeln:

x_next = x +  dt * dx;

Ach verdammt, dann waere man nicht mehr containerunabhaengig. Auch scheint dieses Beispiel nicht mit std::list als Container zu funktionieren, da ihr random access auf die Elemente benoetigt ... getestet habe ich es nicht.

Der grosste Nachteil ist, daß man die Anzahl der Elemente zwischen zwei Iteratoren nicht notwendigerweise mit O(0) rauskriegt.

Wenn man die Methode size einer Liste aufruft, dann garantiert der Standard auch keine O(1) Komplexitaet, sondern auch nur O(n). Moderne Listenimplementationen von std::list achten aber auf ihre Groesse, so dass es zumindest fuer diesen Datentyp keine Probleme gibt. Bei RandomAccessIteratorn sollte es moeglich sein, die Anzahl der Elemente dazwischen in O(1) zu bestimmen.

Das Klassendesign mit virtuellen Funktionen ist Faktor 2 langsamer

Das kann gut sein, wenn es sich um eher "kleine" Differentialgleichungen handelt. Eine virtuelle Funktion in Einfachvererbung kostet halt eine Zeigerdereferenzierung mehr. Ist das im Vergleich zum Rest der Funktion viel, so kann sich das nachteilig auf die Gesamtperformance niederschlagen. Die Funktion wird ja staendig ausgewertet. Die Kosten sind aber pro Funktionsauswertung konstant. Auch hat der Compiler beim Templateeinsatz mehr Moeglichkeiten zu optimieren.

Ich werde mal etwas mit dem neuen Ansatz experimentieren.

headmyshoulder

knivil schrieb:

Sieht interessant aus, vom Ansatz her. Ausprobiert habe ich den Code noch nicht, nur gelesen. Aber etwas zum Codestyle: Beginnt man nach einer geschweiften Klammer eine neue Zeile, so rueckt man ein. Das macht das Lesen angenehmer. Ihr koennt euch ja mal einige Syleguides durchlesen.

Die aktuellen Sourcen befinden sich in https://boost.org/svn/boost/sandbox/odeint nicht daß es da Verwechslungen mit dem alten Code gibt. Dort haben wir auch auf vernünftige Konvertierung geachtet, hoffe ich zumindest.

knivil schrieb:

Auch solltet ihr euch um geeignete Matrix und Vektorklassen bemuehen, die solche Zeilen kapseln:
for( size_t i=0 ; i<dxdt.size() ; ++i ) state[i] += dt * dxdt[i];
und in etwas Vertrauteres verwandeln:
x_next = x +  dt * dx;
Ach verdammt, dann waere man nicht mehr containerunabhaengig. Auch scheint dieses Beispiel nicht mit std::list als Container zu funktionieren, da ihr random access auf die Elemente benoetigt ... getestet habe ich es nicht.

Der code in der boost sandbox sollte mit std::list laufen, man kann halt in der ODE Funktion (oder Klasse) dann keine RandomAccess-Iteratoren benutzen.

knivil schrieb:

Wenn man die Methode size einer Liste aufruft, dann garantiert der Standard auch keine O(1) Komplexitaet, sondern auch nur O(n). Moderne Listenimplementationen von std::list achten aber auf ihre Groesse, so dass es zumindest fuer diesen Datentyp keine Probleme gibt. Bei RandomAccessIteratorn sollte es moeglich sein, die Anzahl der Elemente dazwischen in O(1) zu bestimmen.

Joa, da hast Du Recht. Ich persönlich würde es am besten finden wenn beide Varianten unterstüzt werden. Also einmal, daß die ODE Klasse mit Containern aufgerufen wird, und einmal mit Iteratoren. Allerdings muss man dann wahrscheinlich von jedem Algorithmus zwei Varianten schreiben, einem der die ODE mit Containern aufruft und einem der die ODE mit iteratoren aufruft, oder man abstrahiert das irgendwie weg.

knivil schrieb:

Ich werde mal etwas mit dem neuen Ansatz experimentieren.

Schön, wenn Du Ergebnisse oder Anregungen/Kritik/Vorschläge hast, lass es uns wissen.

knivil

Jap, hatte die anderen Sourcen aus dem uebergeordneten Ordner mit den Speedtests.

PS: Ich habe gerade festgestellt, dass die Welt klein ist. Ich bin auch oefters im Goldenen Kaefig von Golm ... Aber ich dachte, dass es bereits ein Programmpaket fuer Statistik, Odes und den ganzen Kram in der Nichtlinearen Dynamik entwickelt wurde.

headmyshoulder

knivil schrieb:

Jap, hatte die anderen Sourcen aus dem uebergeordneten Ordner mit den Speedtests.

PS: Ich habe gerade festgestellt, dass die Welt klein ist. Ich bin auch oefters im Goldenen Kaefig von Golm ... Aber ich dachte, dass es bereits ein Programmpaket fuer Statistik, Odes und den ganzen Kram in der Nichtlinearen Dynamik entwickelt wurde.

Hehe, ja manchmal ist die Welt wirklich klein.

Ich glaub, hier bei uns wurden schon mehrere solcher und ähnlicher Pakete entwickelt. In der boost gibt es sowas noch nicht, aber ich denke, die Chancen das dort reinzubekommen sind nicht so schlecht. Es wurde mal eine Zeitreihenbibliothek entwickelt, aber die hat es nicht bis zur Veröffentlichung geschafft. Naja, der Sinn ist auch endlich mal was zu schaffen, was Standards benutzt. Da gibt es echt nicht viel. Vielleicht benutzen dann auch wieder mehr Leute C++ statt Python. Das verbreitet sich in der Wissenschaft gerade rasant^^.

FelixKlein

Hallo headmyshoulder!

Ich habe erst gestern das entsprechende Kapitel in den Num. Recipes gelesen und bin heute auf deine Idee von dem C++ Odeint Framework gestoßen. Ich hoffe das Projekt läuft noch!

Ich selber habe Mathematik studiert. Mit der ODE-Solvern und Schrittweitensteurungen beschäftige mich erst seit Kurzem. Auch in C++ bin ich noch relativ unerfahren. Aber ich werde muss mich sowieso auf beiden Gebieten in nächster Zeit einarbeiten. Über jeden Austausch und gegebenenfalls eine Zusammenarbeit würde ich mich freuen!

Viele Grüße
FelixKlein

headmyshoulder

Hi FelixKlein,

das odeint Projekt läuft noch und ist auch schon recht fortgeschritten, du hast quasi eine menge aufzuholen. Zu finden ist das Ganze in der boost sandbox

https://www.boost.org/svn/boost/sandbox/odeint

Akutell könnten wir die Dormand-Prince Methoden gebrauchen. Wenn Du Bock kannst Du dir ja ansehen wie unsere lib funktioneirt und dann die Methoden implementieren.

Viele Grüße,

headmyshoulder