[R] Numerik und Template-Metaprogrammierung

headmyshoulder

So, ich hab den Artikel jetzt auf R gesetzt. Ich kann natürlich trotzdem "fachliche" Anmerkungen einfügen. Ich wollte auch die Sourcen hochladen, aber irgednwie krieg ich keine Verbindung auf den ftp Server. Vielleicht sind Username/PW nicht mehr aktuell oder der ftp ist gerade down.

GPC

headmyshoulder schrieb:

So, ich hab den Artikel jetzt auf R gesetzt. Ich kann natürlich trotzdem "fachliche" Anmerkungen einfügen. Ich wollte auch die Sourcen hochladen, aber irgednwie krieg ich keine Verbindung auf den ftp Server. Vielleicht sind Username/PW nicht mehr aktuell oder der ftp ist gerade down.

Also ich habe gerade eben erfolgreich eine Verbindung auf den server hinbekommen. Vermutlich hast du dich vertippt

Und den Artikel werd ich heute oder morgen noch korrekturlesen. Denk auch dran, dass du noch ein Autorenprofil brauchst, in dem du ein wenig über dich schreibst. Muss ja nix großes sein

headmyshoulder

GPC schrieb:

Also ich habe gerade eben erfolgreich eine Verbindung auf den server hinbekommen. Vermutlich hast du dich vertippt

Nimmst Du die Daten von hier? Und wenn ja, normales FTP ohne verschlüsselung? Und welcher Port? Vielleicht ist mein Client auch nicht so intelligent (filezilla).

GPC schrieb:

Und den Artikel werd ich heute oder morgen noch korrekturlesen. Denk auch dran, dass du noch ein Autorenprofil brauchst, in dem du ein wenig über dich schreibst. Muss ja nix großes sein

Jepp, mach ich.

GPC

headmyshoulder schrieb:

GPC schrieb:

Also ich habe gerade eben erfolgreich eine Verbindung auf den server hinbekommen. Vermutlich hast du dich vertippt

Nimmst Du die Daten von hier? Und wenn ja, normales FTP ohne verschlüsselung? Und welcher Port? Vielleicht ist mein Client auch nicht so intelligent (filezilla).

Jappa, nehm ich. Und standardport. Hab's mit dem ftp-Konsolenprogramm gemacht

gpc@desktop:~$ ftp c-plusplus.net
Connected to c-plusplus.net.
220 (vsFTPd 2.0.7)
Name (c-plusplus.net:gpc): magazin_public
331 Please specify the password.
Password:
230 Login successful.
Remote system type is UNIX.
Using binary mode to transfer files.
ftp> ls
200 PORT command successful. Consider using PASV.
150 Here comes the directory listing.
drwxr-sr-x    2 1003     33           4096 Dec 21  2007 Compiler
drwxrwsr-x    2 33       33           4096 Oct 03  2005 Debug_VC6
[.....]

Kann's sein dass du ein O (Buchstabe) anstatt der Ziffer 0 nimmst (direkt nach pE1)?

headmyshoulder

GPC schrieb:

Jappa, nehm ich. Und standardport. Hab's mit dem ftp-Konsolenprogramm gemacht

Danke, mit dem Konsolenprogramm komm ich jetzt rein. Allerdings hab ich jetzt das Problem das ich kein Verzeichnis anlegen kann:

mkdir template_metaprogrammierung_numerik

liefert

550 Create directory operation failed.

Ansonsten seh ich alle Verzeichnisse und kann mir auch den Inhalt ansehen. Sry, falls ich mich zu doof anstelle.

GPC

headmyshoulder schrieb:

GPC schrieb:

Jappa, nehm ich. Und standardport. Hab's mit dem ftp-Konsolenprogramm gemacht

Danke, mit dem Konsolenprogramm komm ich jetzt rein. Allerdings hab ich jetzt das Problem das ich kein Verzeichnis anlegen kann:
mkdir template_metaprogrammierung_numerik
liefert
550 Create directory operation failed.
Ansonsten seh ich alle Verzeichnisse und kann mir auch den Inhalt ansehen. Sry, falls ich mich zu doof anstelle.

Sorry für die späte Antwort, die Woche war etwas stressig. Man kann atm tatsächlich kein Unterverzeichnis anlegen, habe da aber intern mal nachgefragt und sobald das wieder geht, geb ich dir hier Bescheid

GPC

Okay, also die Berechtigungen sind jetzt frei. Ich hab das Verzeichnis schon mal unter dem Namen "template_metaprogrammierung_numerik" angelegt, du kannst jetzt die Inhalte selbst hochladen

GPC

Einführung

In diesem Artikel werde ich einige Techniken im Zusammenhang mit Templatemetaprogrammierung und numerischen Verfahren vorstellen. Insbesondere wird es dabei um numerische Verfahren zum Lösen von Differentialgleichungen (DGL) gehen. Die hier vorgestellten Methoden lassen sich aber problemlos auf eine Vielzahl von anderen numerischen Problemen anwenden.

Templatemetaprogrammierung ist eine moderne Programmiertechnik bei welcher der C++ Compiler elementare Berechnungen durchführt oder komplexen C++ Code und Algorithmen erzeugt. Ein entscheidender Aspekt dabei ist, dass Templatemetaprogramme während des Compiliervorganges ausgewertet bzw. aufgeführt werden. Das klassische Hello world-Programm ist hier die Berechnung der Fakultät, welche direkt vom Compiler berechnet wird. Zur Laufzeit des fertigen Programmes wird diese Berechnung keine Zeit in Anspruch nehmen.

Templatemetaprogrammierung kann auch für komplexe Probleme und Algorithmen verwendet werden. Ein Beispiel ist die Bibliothek Boost.MetaStateMachine, mit welcher Zustandsmachinen anhand einer (statischen) Übergangstabelle während des Kompilierens erzeugt werden können. Ein weiteres Beispiel ist Boost.Accumulators, eine Bibliothek für statistische Auswertungen. Dort werden die genauen Berechnung und statistischen Routinen auch während der Kompilierphase erzeugt. Ausserdem werden Techniken der Templatemetaprogrammierung auch oft im Zusammenhang mit Typkonvertierungen und Typüberprüfungen benutzt.

Es gibt eine Vielzahl von Bibliotheken welche mehr oder weniger mit dem Thema Templatemetaprgrammierung zu tun haben. Die meisten davon sind in den Boost Bibliotheken zu finden. Bevor man anfängt Metaprogramme selbst zu schreiben und zu implementieren, empfehle ich aber zuerst dort nachzusehen, ob nicht bereits eine adäquate Lösung existiert. Selbst durch das konsequente Anwenden und Benutzen dieser Bibliotheken kann man sehr viel lernen. Als wichtigste Bibliothek ist dabei sicherlich Boost.MPL mit einer Vielzahl von Metafunktionen- und algorithmen und Compile-time Sequenzen. Eine weitere wichtige Bibliothek ist Boost.Fusion, in welcher die Runtime-Aspekte mehr betont werden. Zusätzlich müssen in der Templatemetaprogrammierung sehr oft Typen klassifiziert werden. Dafür stehen etliche Funktion in Boost.TypeTraits zur Verfügung. Beispielsweise werden Metafunktionen wie is_pointer , is_class , aber auch Konversionen add_const , remove_const zur Verfügung gestellt. Ein Teil dieser Funktionalität wird auch Einzug in den kommenden C++ Kompilerstandard finden.

Zur Entwicklung von Meta-Algorithmen werden wir hier ausschliesslich Boost.MPL verwendet. Diese Bibliothek ist in ihrem Umfang und Anwendungsmöglichkeit am Weitesten entwickelt. Es gibt schlicht und einfach keine wirklichen Alternativen dazu.

Wer sich ernsthaft mit dem Thema Templatemetaprogrammierung auseinandersetzen will, dem würde ich einige Pflichtlektüren empfehlen. "Modern C++ Design" von Andrei Alexandrescu ist vielleicht eine guter, wenn auch nicht einfacher Einstieg. Dort werden sogut wie alle der grundlegenden Techniken erklärt, aber man sollte schon recht sicher im Umgang mit C++ sein. Das Buch "C++ Template Metaprogramming" von David Abrahams und Aleksey Gurtovoy erklärt die Thematik recht anschaulich am Beispiel von Boost.MPL, ohne allerdings detailiert die Implementierungen der einzelnen Algorithemen zu erläutern. Und zu guter Letzt sei noch das Buch "C++ Templates" von David Vandevoorde und Nicolai Josuttis erwähnt, welches sich ausführlichst mit Templates beschäftigt und damit die (theoretischen) Grundlagen für vielen Methoden liefert.

Der Artikel ist wie folgt gegliedert. Nach dieser Einführung werde ich im zweiten Teil zeigen, wie numerische Konstanten für Templatemetaprogrammierung definiert werden. Wir werden diese Definition später, wenn es um die Konstruktion der numerische Algorithmen geht, wieder benötigen. Danach werde ich anhand eines Beispieles kurz erläutern, welche Probleme generell beim Erstellen von Templatemetaprogrammen auftreten können, wie man am besten Fehler in diesen Programmen findet und Templatemetaprogramme debuggt. Danach kommt der Hauptteil, in welchem die klassischen Runge-Kutta-Verfahren zum Lösen von Differentialgleichungen mit Hilfe von Templatemetaprogrammen konstruiert werden. Der Artikel endet mit einer Diskussion der Ergebnisse und Methoden.

Numerische Konstanten

Numerische Konstanten können mit den eingebauten integralen Typen der MPL dargestellt werden:

typedef mpl::int_< 0 > zero;
typedef mpl::char_< 1 > one;
typedef mpl::long_< 2 > two;
typedef mpl::long_< -2 > minus_two;

Der numerische Wert wird gewöhnlich in einer statischen Membervariable ::value abgelegt und kann wie folgt ausgegeben werden.

cout << zero::value << endl;
cout << one::value << endl;
cout << two::value << endl;
cout << minus_two::value << endl;

Ausserdem sind diese Typen nullary Metafunktionen, welche ihren Typ zurückgeben. Daher sind die beiden Typen mpl::int_< 1 > und mpl::int_< 1 >::type gleichwertig.

Man beachte: Da Zahlen als Typen definiert sind, steht ihr Wert während des Kompiliervorgangs fest. Ausdrücke, welche diese Werte benutzen, können dann vom Kompiler optimiert werden. Allerdings können nur die ganzzahligen Typen (char,short,int,long,...) mit Hilfe der numerischen Konstanten der mpl ausgedrückt werden. Benötigen wir Fliesskommazahlen zur Compilezeit, so können wir diese durch rationale Zahlen darstellen:

template< class T , T Numerator , T Denominator >
struct rational_number
{
    typedef T value_type;
    typedef rational_number type;
    const static value_type numerator = Numerator;
    const static value_type denominator = Denominator;
};

template< class Num >
struct conv
{
    static double get_value( void )
    {
        return static_cast< double >( Num::value );
    }
};

template< class T , T Numerator , T Denominator >
struct conv< rational_number< T , Numerator , Denominator > >
{
    typedef rational_number< T , Numerator , Denominator > num_type;
    static double get_value( void )
    {
        return static_cast< double >( num_type::numerator ) / static_cast< double >( num_type::denominator );
    }
};

Die Klasse conv übernimmt hier nur die Konvertierung zu double s. Man sieht aber hier auch einen bedeutenden Unterschied zu integralen Typen. Der Wert einer rationalen Zahl kann nicht so ohne Weiteres in einer statischen Member-Variable abgelegt, werden, da dies nur für integrale Typen erlaubt ist. Es sein denn, man definiert diese Variable extern, was aber nicht besonders ansprechend ist. Ausserdem ist die Variable dann auch wirklich als Symbol vorhanden. Daher muss hier der Umweg über eine statische Klasse gegangen werden. Diese Klasse ist spezialisiert für rational_number und man kann mit ihr sowohl integrale als auch Fließkommazahlen einheitlich benutzen, beispielsweise so:

typedef rational_number< long , 10 , 15 > two_third;
cout << conv< two_third >::get_value() << endl;
cout << conv< two >::get_value() << endl;

Für unsere Anwendung ist die Definition von numerischen Konstanten völlig ausreichend. Zusätzlich könnte man jetzt noch die Operation (+,-,*,/) definieren. Es gibt sogar einige Bibliotheken, welche die Benutzung von numerischen Konstanten soweit auf die Spitze treiben, dass beispielsweise trigonometrische Funktionen wie der Sinus oder der Cosinus vom Compiler berechnet werden können, siehe hier.

Beispielalgorithmen und Debugging

Debugging von Templatemetaalgorithmen ist meistens problematisch. Eine richtige Lösung hab ich nicht, allerdings kann man mit Hilfe von typeid einige Informationen über mögliche Probleme erhalten. Nehmen wir ein Beispiel: Wir wollen einen beliebig langen fusion Vektor konstruieren, dessen Elemente tr1::array s sind. Die Länge der Arrays soll dabei ihrem Index in dem Vektor entsprechen. Beispielsweise soll ein Vektor der Länge 3 so aussehen

typedef fusion::vector
<
    tr1::array< double , 1 > ,
    tr1::array< double , 2 > ,
    tr1::array< double , 3 >
> coef_type;

Ein erster Versuch solch einen Vektor zu konstruieren wäre

typedef mpl::copy
<
    mpl::range_c< size_t , 1 , N + 1 > ,
    mpl::inserter
    <
        mpl::vector<> , 
        mpl::push_back
        <
            mpl::_1 ,
            tr1::array< double , mpl::_2::value >
        >
    >
>::type mpl_coef_type;
typedef fusion::result_of::as_vector< mpl_coef_type >::type coef_type;

Dies würde einem klassischen C++ Kopieralgorithmus mit einem speziellen back_insert_iterator entsprechen, ganz analog zu

// Achtung, das Beispiel wird nicht kompilieren
vector< size_t > indices;
for( size_t i=1 ; i<=N ; ++i ) indices.push_back( i );
vector< arrays > vec;
copy( indices.begin() , indices.begin() , back_insert_iterator< vector< arrays > >( vec ) );

Aber gehen wir Schritt für Schritt durch den Meta-Algorithmus:

mpl::copy< Seq , Inserter > ist der zentrale Kopiervorgang. Er nimmt eine MPL-Sequenz Seq entgegen und gibt jedes Element dieser Sequenz an den Inserter weiter.
mpl::range_c< size_t , 1 , N + 1 > konstruiert einen mpl::vector dessen Elemente size_t s von 1 bis N sind, also wie mpl::vector< 1 , 2 , 3 , ... , N-1 , N > . Dieser Vektor wird benötigt, da er erstens die richtige Anzahl von Elementen für den finalen Vektor besitzt und zweitens jedes Element die Länge des aktuellen Arrays ist.
mpl::inserter< Seq, Op > ist ein allgemeiner Inserter, ähnlich zu den Insert-Iteratoren der STL. Als erstes Template-Argument wird hier eine initiale Sequenz erwartet, in unserem Falle ein leerer Vektor mpl::vector<> . Das zweite Argument definiert die Einfüge-Operation.
mpl::push_back< Seq , El > macht nichts anderes als El an die Sequenz Seq anzuhängen. In unserem Fall wird also an den aktuellen Vektor ein neues Element vom Typ std::tr1::array angehängt. Hier sieht man auch sehr schön die Verwendung von Platzhaltern. mpl::_1 und mpl::_2 bezeichnen das erste bzw. zweite Argument, welches an push_back übergeben wird.
fusion::result_of::as_vector konvertiert den Mpl-Vektor in einen Fusion-Vektor.

Ausgeschrieben soll wird der gesamte Algorithmus in etwa so aussehen:

mpl::push_back
<
    mpl::push_back
    < 
        mpl::push_back
        <
            mpl::vector<> ,
            tr1::array< double , 1 >
        >::type ,
        tr1::array< double , 2 >
    >::type ,
    tr1::array< double , 3 >
>::type

Der Code kompiliert und man könnte meinen, dass wir die Aufgabe geschafft haben. Aber dem ist nicht so. Wollen wir diesen jetzt verwenden, wird der Compiler eine schier unglaubliche Anzahl von Fehlern ausspucken. Diese nach ihrer Essenz zu durchforsten ist mühsam und selbst wenn man die entscheidenden Stellen gefunden hat, kann es sein dass man mit der Fehlermeldung nicht viel anfangen kann (mir geht das zumindest so). Um diese Fehlersuche abzukürzen gibt es mehrere Wege:

Benutze STLFilt! Es wird die Anzahl der Fehlermeldungen drastisch reduzieren.
Benutze statische Assertions! Damit können bestimmte Fehler viel früher gefunden werden. Statische Assertions sind insbesondere für Library-Entwickler wichtig, weil damit Bedingungen an die Typen geprüft werden.
Benutze einen Typdebugger! Er wird Information über den Typ liefern, welche eventuell hilfreich sind.

template< class T >
struct wrap {};

struct debug_type
{
    size_t m_indent;

    debug_type( size_t indent = 0 )
    : m_indent( indent ) { }

    template< class T >
    void operator()( wrap<T> ) const
    {
        typedef typename boost::mpl::is_sequence< T >::type is_seq;
        debug_impl( wrap< T >() , is_seq() );
    }

    template< class Seq >
    void debug_impl( wrap< Seq > , boost::mpl::true_ ) const
    {
        std::cout << indent() << "Sequence : " << typeid(Seq).name() << std::endl;
        boost::mpl::for_each<
        typename boost::mpl::transform< Seq , wrap< boost::mpl::_1 > >::type >( debug_type( m_indent + 1 ) );
    }

    template< class T >
    void debug_impl( wrap< T > , boost::mpl::false_ ) const
    {
        std::cout << indent() << typeid(T).name() << std::endl;
    }

    std::string indent( void ) const
    {
        std::string in;
        for( size_t i=0 ; i<m_indent ; ++i ) in += "  ";
        return in;
    }
};

template< class T >
void debug( size_t indent = 0 )
{
    debug_type d( indent );
    wrap< T > t;
    d( t );
}

Der wrap -Trick kommt aus dem Buch von David Abrahams und sorgt dafür, dass der Typ-Printer auch mit Referenzen funktioniert (Referenzen sind nicht default-konstruierbar und können nicht einfach von mpl::for_each erzeugt werden).

In unserem Beispiel des Vektors von Arrays führt nun die Benutzung von

debug< mpl_coef_type >();
debug< coef_type >();
dout << endl;

zu folgender Ausgabe:

Sequence : N5boost3mpl6v_itemINSt3tr15arrayIdLm2EEENS1_IS4_NS1_IS4_NS0_6vectorIN4mpl_2naES7_S7_S7_S7_S7_S7_S7_S7_S7_S7_S7_S7_S7_S7_S7_S7_S7_S7_S7_EELi0EEELi0EEELi0EEE
  NSt3tr15arrayIdLm2EEE
  NSt3tr15arrayIdLm2EEE
  NSt3tr15arrayIdLm2EEE
Sequence : N5boost6fusion7vector3INSt3tr15arrayIdLm2EEES4_S4_EE
  NSt3tr15arrayIdLm2EEE
  NSt3tr15arrayIdLm2EEE
  NSt3tr15arrayIdLm2EEE

Und Voila, man sieht dass mpl::_2::value nicht richtig ausgelöst/ersetzt wurde, denn die Arrays haben die Länge 2. Eine Lösung ist schnell gefunden, wir packen die Definition des Arrays einfach in eine Meta-Funktion:

template< class T , class Constant >
struct array_wrapper
{
    typedef typename tr1::array< T , Constant::value > type;
};

typedef mpl::copy
<
    mpl::range_c< size_t , 1 , N + 1 > ,
    mpl::inserter
    <
        mpl::vector<> ,
        mpl::push_back
        <
            mpl::_1 ,
            array_wrapper< double , mpl::_2 >
        >
    >
>::type mpl_coef_type;
typedef fusion::result_of::as_vector< mpl_coef_type >::type coef_type;

und unser kleiner aber feiner Algorithmus liefert das gewünschte Ergebnis.

Konstruktion von numerischen Algorithmen

Kommen wir nun zum eigentlich Thema dieses Artikels, der Konstruktion von numerischen Algorithmen mit Hilfe von Template-Metaprogramming. Wir verdeutlichen unseren Ansatz anhand der klassischen expliziten Runge-Kutta-Solvern zum Lösen von gewöhnlichen Differentialgleichungen. Ähnliche Ansätze können auch für verwandte numerische Probleme wie

Numerische Integration / Lagrange Integration,
Konstruktion von Spline-Basis-Funktion, oder auch
Finite Differenzen zum Lösen von partiellen Differentialgleichungen

genutzt werden.

Ein gewöhnliche Differentialgleichung beschreibt die (zeitliche) Entwicklung einer Größe x anhand ihrer Ableitung(en):

dx / dt = f(x,t)

x kann dabei vektoriell sein, daher aus mehreren Werten bestehen. Existiert zu einem Zeitpunkt t=t0 eine Anfangsbedingung x(t0)=x0 und ist die zeitliche Entwicklung von x gesucht, so spricht man von einem Anfangswertproblem. Numerisch werden Anfangswertprobleme normalerweise iterativ gelöst. Man bestimmt also eine aufeinanderfolgende Sequenz x(t0),x(t0+dt),x(t0+2dt),... welche die Lösung numerisch repäsentiert. dt ist dabei die Schrittweite.

Eine einfache Klasse von solchen Algorithmen sind sogenannte explizite Runge-Kutta-Verfahren. Sie haben den grossen Vorteil, dass keine Gleichungssysteme gelöst werden müssen. Schematisch werden Runge-Kutta-Solver durch Butcher-Tableaus dargestellt:

c[t]1[/t] |
c[t]2[/t] | a[t]2,1[/t]
c[t]3[/t] | a[t]3,1[/t] a[t]3,2[/t]
...
c[t]s[/t] | a[t]s,1[/t] a[t]s,2[/t] ... a[t]s,s-1[/t]
-----------------------------------
   | b[t]1[/t] b[t]2[/t] ... b[t]s-1[/t] b[t]s[/t]

s bezeichnet dabei die Ordnung des Verfahrens, was ein wichtiger Parameter ist, da er die Genauigkeit der Lösung charakterisiert.

Aus solch einer Tabelle kann man nun das Verfahren zum Lösen der DGL wie folgt konstruieren

k[t]1[/t] = f( x(t) , t + c[t]1[/t] * dt )
k[t]2[/t] = f( x(t) + dt * a[t]2,1[/t] * k[t]1[/t] , t + c[t]2[/t] * dt )
k[t]3[/t] = f( x(t) + dt * ( a[t]3,1[/t] * k[t]1[/t] + a[t]3,2[/t] * k[t]2[/t] ), t + c[t]3[/t] * dt )
...
k[t]s[/t] = f( x(t) + dt * ( a[t]s,1[/t] * k[t]1[/t] + ... + a[t]s,s-1[/t] * k[t]s-1[/t] ) , t + c[t]s[/t] * dt )
x(t+dt) = x(t) + dt * ( b[t]1[/t] k[t]1[/t] + ... + b[t]s[/t] k[t]s[/t] )

Ein einfaches Beispiel ist das Euler-Verfahren:

0 |
--|---
  | 1

Oder ausgeschrieben:

k[t]1[/t] = f( x(t) , t )
x(t+dt) = x(t) + dt * k[t]1[/t]

Ein etwas komplizierteres Beispiel ist das klassischen Runge-Kutta-Verfahren vierter Ordnung (im folgenden als RK4 bezeichnet):

0   |
1/2 | 1/2
1/2 | 0   1/2
1   | 0   0   1
----|----------------
    | 1/3 1/6 1/6 1/3

In Pseudocode sieht diese Methode so aus:

k[t]1[/t] = f( x(t) , t )
k[t]2[/t] = f( x(t) + dt * 1/2 * k[t]1[/t] , t + 1/2 * dt )
k[t]3[/t] = f( x(t) + dt * 1/2 * k[t]2[/t] , t + 1/2 * dt )
k[t]4[/t] = f( x(t) + dt * 1 * k[t]3[/t] , t + 1 * dt )
x(t+dt) = x(t) + dt * ( 1/6 k[t]1[/t] + 1/3 k[t]2[/t] + 1/3 k[t]3[/t] + 1/6 k[t]4[/t] )

Wie werden nun Schritt für Schritt einen Meta-Algorithmus konstruieren, welcher solche Algorithmen während des Kompiliervorganges erzeugt. Als erstes werden die Koeffizienten a[t]i,j[/t] , b[t]i[/t] und c[t]i[/t] definiert. Danach wird das Hauptgerüst des Algorithmus entworfen und im letzten Teil die einzelnen Schritte ausgeführt. Zum Schluss wird ein alternative Implementation vorgestellt, welche sich durch eine einfache Konstrukion der Koeffizienten auszeichnet.

Die Koeffizienten eines Verfahrens werden in mpl::vector abgelegt. Für RK4 sieht die Definition wie folgt aus:

// allgemeine Konstanten
typedef mpl::int_< 0 > null;
typedef mpl::int_< 1 > one;
typedef rational_number< long , 1 , 2 > one_half;
typedef rational_number< long , 1 , 3 > one_third;
typedef rational_number< long , 1 , 6 > one_sixth;

typedef mpl::vector
<
    mpl::vector< one_half > ,
    mpl::vector< null , one_half > ,
    mpl::vector< null , null , one >
> rk4_a;
typedef mpl::vector< null , one_half , one_half , one > rk4_c;
typedef mpl::vector< one_sixth , one_third , one_third , one_sixth > rk4_b;

Der Solver selber ist eine templatisierte Klasse, welcher durch die Koeffizienten parametrisiert ist:

template< class State , class ButcherA , class ButcherB , class ButcherC >
class explicit_runge_kutta
{
public:
    typedef ButcherA butcher_a;
    typedef ButcherB butcher_b;
    typedef ButcherC butcher_c;

    static const size_t num_of_sages = mpl::size< ButcherC >::value;

    typedef double time_type;
    typedef State state_type;
...
}

time_type und state_type geben die Typen von t bzw. x an, und für state_type wird ein std::tr1::array< double , N > erwartet.

Als nächsten überprüfen wir, ob butcher_a , butcher_b and butcher_c den Anforderungen an die Koeffizienten genügen. Insbesondere prüfen wir, ob diese Typlisten die entsprechenden Anzahlen von Elementen enthalten.

template< class State , class ButcherA , class ButcherB , class ButcherC >
class explicit_runge_kutta
{
    ...
    BOOST_STATIC_ASSERT( ( num_of_stages > 0 ) );
    BOOST_STATIC_ASSERT( ( num_of_stages == size_t( mpl::size< butcher_b >::value ) ) );
    BOOST_STATIC_ASSERT( ( ( num_of_stages - 1 ) == size_t( mpl::size< butcher_a >::value ) ) ) ;

    typedef typename mpl::equal
    <
        mpl::range_c< int , 1 , num_of_stages > ,
        typename mpl::transform< butcher_a , mpl::size< mpl::_1 > >::type ,
        mpl::equal_to< mpl::_1 , mpl::_2 >
    >::type equal_type;
    BOOST_STATIC_ASSERT(( equal_type::value == true ));
};

Jetzt kommen wir zum eigentlichen Kernstück unseres Metaalgorithmus. Die expliziten Runge-Kutta-Verfahren verfügen über s Stufen, wobei s hier wieder die Ordnung des Verfahrens ist. Die erste Stufe sieht wie folgt aus

k[t]1[/t] = f( x(t) , t + c[t]1[/t] )
xtmp = x(t) + dt * a[t]2,1[/t] * k[t]1[/t]

In den nachfolgenden Stufen (ausser der letzten ) wird immer

k[t]i[/t] = f( xtmp , t + c[t]i[/t] )
xtmp = x(t) + dt * a[t]i+1,1[/t] * k[t]1[/t] + ... + a[t]i+1,i[/t] * k[t]i[/t]

berechnet. Auch die letzt Stufe unterscheidet sich, hier wird das Ergebnis geschrieben

k[t]s[/t] = f( xtmp , t + c[t]s[/t] )
x(t+dt) = x(t) + dt * b[t]1[/t] * k[t]1[/t] + ... + b[t]s[/t] * k[t]s[/t]

Diesen Algorithmus setzten wir nun in einen Meta-Algorithmus um:

template< class State , class ButcherA , class ButcherB , class ButcherC >
class explicit_runge_kutta_mpl
{
public:
    ...
    typedef mpl::range_c< size_t , 0 , num_of_stages > indices;
    typedef typename mpl::push_back< butcher_a , butcher_b >::type butcher_right;
    typedef mpl::zip_view< mpl::vector< indices , butcher_c , butcher_right > > butcher_tableau;

    template< class System >
    void do_step( System &system , state_type &x , double t , const double dt )
    {
        mpl::for_each< butcher_tableau >( calculate_stage< System >( system , x , m_x_tmp , m_k_vector , t , dt ) );
    }
private:

    state_type m_k_vector[num_of_stages];
    state_type m_x_tmp;
    ...
};

Schauen wir uns an, was hier passiert. butcher_tableau ist ein Mpl-Vektor mit folgenden Einträgen:

0 , c[t]1[/t] , a[t]2,1[/t] 
1 , c[t]2[/t] , a[t]3,1[/t] , a[t]3,2[/t]
2 , c[t]3[/t] , a[t]4,1[/t] , a[t]4,2[/t] , a[t]4,3[/t] 
...
s , c[t]s[/t] , b[t]1[/t]   , b[t]2[/t]   , b[t]3[/t]   , ... , b[t]s[/t]

Die i.te Zeile enthält also alle Informationen, welche wir zur Berechnung der i.ten Stufe benötigen. Über diesen Vektor wird nun in der Zeile

mpl::for_each< butcher_tableau >( calculate_stage< System >( system , x , m_x_tmp , m_k_vector , t , dt ) );

iteriert und in jeder Iteration calculate_state aufgerufen. Die Methode do_step geht nun einen Zeitschritt in der Differentialgleichung nach vorn, do_step berechnet x(t+dt) anhand von x(t) . Die Differentialgleichung selber wird durch system angegeben. system kann hier ein Funktor aber auch ein Funktionspointer sein. Es wird dabei immer erwartet, dass wir

state_type dxdt;
system( x , dxdt , t );

ausführen können.

Nun wenden wir uns der der eigentlichen Berechnung jeder einzelnen Stufe zu. Das passiert in calculate_stage , welche eine Member-Struktur von explicit_runge_kutta_mpl ist. Diese Struktur ist gleichzeitig auch ein Funktor, daher existiert der Klammer-Operator.

template< class State , class ButcherA , class ButcherB , class ButcherC >
class explicit_runge_kutta_mpl
{
public:
    ...
    template< class System >
    struct calculate_stage
    {
        System &system;
        state_type &x , &x_tmp;
        state_type *k_vector;
        const double t;
        const double dt;

        calculate_stage( System &_system , state_type &_x , state_type &_x_tmp , state_type *_k_vector , const double _t , const double _dt )
        : system( _system ) , x( _x ) , x_tmp( _x_tmp ) , k_vector( _k_vector ) , t( _t ) , dt( _dt )
        {}

        template< class Stage >
        void operator()( Stage v )
        {
            typedef typename mpl::at< Stage , mpl::int_< 0 > >::type index_type;
            typedef typename mpl::at< Stage , mpl::int_< 1 > >::type time_value_type;
            typedef typename mpl::at< Stage , mpl::int_< 2 > >::type coef_type;

            const static size_t index = index_type::value;
            const static double time_value = conv< time_value_type >::get_value();

            BOOST_STATIC_ASSERT(( ( mpl::size< coef_type >::value - 1 ) == index ));

            if( index == 0 ) system( x , k_vector[index] , t + time_value * dt );
            else system( x_tmp , k_vector[index] , t + time_value * dt );

            if( index == ( num_of_stages - 1 ) )
                algebra< state_type , coef_type , mpl::int_< index > >::do_step( x , x , k_vector , dt );
            else
                algebra< state_type , coef_type , mpl::int_< index > >::do_step( x_tmp , x , k_vector , dt );
        }
    };
    ...
};

Das Klassen-Template algebra enthält nur eine einzelne statische Methode do_step und ist durch den Zustandstyp state_type , die Koeffizienten des Verfahren coef_type und die Ordnung des Verfahrens order parametrisiert.

template< class state_type , class coef_tye , class order > struct algebra_mpl;

algebra muss für alle Ordnungen spezialisiert werden. Hier könnte bestimmt der Präprozessor helfen, aber das ist eine andere Geschichte. Deswegen schreiben wir alle Ordnungen aus. Zum Beispiel sehen die ersten beiden Ordnungen so aus:

template< class state_type , class coef_type >
struct algebra_mpl< state_type , coef_type , mpl::int_< 0 > >
{
    static void do_step( state_type &x_tmp , const state_type &x , const state_type *k_vector , const double dt )
    {
        const static double a1 = dt * conv< typename mpl::at< coef_type , mpl::int_< 0 > >::type >::get_value();

        for( size_t i=0 ; i<x.size() ; ++i )
        {
            x_tmp[i] = x[i] + a1 * k_vector[0][i];
        }
    }
};

template< class state_type , class coef_type >
struct algebra_mpl< state_type , coef_type , mpl::int_< 1 > >
{
    static void do_step( state_type &x_tmp , const state_type &x , const state_type *k_vector , const double dt )
    {
        const static double a1 = dt * conv< typename mpl::at< coef_type , mpl::int_< 0 > >::type >::get_value();
        const static double a2 = dt * conv< typename mpl::at< coef_type , mpl::int_< 1 > >::type >::get_value();

        for( size_t i=0 ; i<x.size() ; ++i )
        {
            x_tmp[i] = x[i] + a1 * k_vector[0][i] + a2 * k_vector[1][i];
        }
    }
};

Damit sind wir fertig. Das letzte was uns jetzt noch übrig bleibt, ist einen speziellen Solver zu initialisieren sowie zu instanzieren und dann können wir mit dem Lösen von Differentialgleichungen beginnen. Für das vorangegangene Euler-Beispiel sieht die Definition des Solver wie folgt aus:

typedef mpl::int_< 0 > null;
typedef mpl::int_< 1 > one;

typedef mpl::vector< > euler_a;
typedef mpl::vector< one > euler_b;
typedef mpl::vector< null > euler_c;

template< class State > class stepper_euler
: public explicit_runge_kutta_mpl< State , euler_a , euler_b , euler_c > {};

Als erstes werden die Konstanten für das Verfahren definiert, in diesem Fall 0 und 1. Danach werden die Koeffizienten a[t]i,j[/t] , b[t]i[/t] und c[t]i[/t] definiert. Die Koeffizienten b und c enthalten jeweils nur einen Eintrag, die Koeffizienten a sind sogar leer. Danach wird der Euler-Solver per Vererbung definiert. Für den klassischen Runge-Kutta-Solver sieht alles sehr ähnlich aus, nur dass andere Koeffizienten vorkommen:

typedef rational_number< int , 1 , 2 > one_half;
typedef rational_number< int , 1 , 3 > one_third;
typedef rational_number< int , 1 , 6 > one_sixth;

typedef mpl::vector
<
    mpl::vector< one_half > ,
    mpl::vector< null , one_half > ,
    mpl::vector< null , null , one >
> rk4_a;
typedef mpl::vector< null , one_half , one_half , one > rk4_c;
typedef mpl::vector< one_sixth , one_third , one_third , one_sixth > rk4_b;

template< class State > class stepper_rk4
: public explicit_runge_kutta_mpl< State , rk4_a , rk4_b , rk4_c > { };

Alternative Implementation

Ein Nachteil der obigen Implementierung des Compile-Zeit Runge-Kutta Algorithmus ist sicherlich die Art und Weise wie die Koeffizienten angegeben werden müssen. Für einen Koeffizienten wie 1/6 müssen wir beispielsweise rational_number< int , 1 , 6 > schreiben. Eine einfachere Definition könnte beispielsweise konstante Arrays benutzen, welche im Konstruktor des Runge-Kutta-Verfahrens übergeben werden. Wir werden diesen Ansatz jetzt entwickeln, er erfordert allerdings ein kleines bisschen mehr Metaprogrammierungszauberei.

In der alternativen Implementierung werden die Typen der Koeffizienten folgendermaßen definiert:

template< class T , class Constant >
struct array_wrapper
{
    typedef typename std::tr1::array< T , Constant::value > type;
};

template< class StateType , size_t NumOfStages >
class explicit_runge_kutta_fusion
{
public:
    typedef StateType state_type;
    const static size_t num_of_stages = NumOfStages;

    typedef typename fusion::result_of::as_vector
    <
        typename mpl::copy
        <
            mpl::range_c< size_t , 1 , num_of_stages > ,
            mpl::inserter
            <
                mpl::vector0< > ,
                mpl::push_back< mpl::_1 , array_wrapper< double , mpl::_2 > >
            >
        >::type
    >::type coef_a_type;

    typedef std::tr1::array< double , num_of_stages > coef_b_type;
    typedef std::tr1::array< double , num_of_stages > coef_c_type;
    ...
}

Die Koeffizienten b und c werden in einfachen Arrays gespeichert. Für die Koeffizienten a benötigen wir allerdings etwas Komplizierteres: einen fusion::vector . Solche Vektoren sind ähnlich zu std::tr1::tuple , aber mit dem Unterschied dass dazu eine riesige Menge an Meta-Algorithmen und Funktionen zur Verfügung stehen. Der Meta-Algorithmus hier konstruiert einen solchen Vektor (siehe auch das obige Beispiel). result_of::as_vector konvertiert einen MPL-Vektor in einen Fusion-Vektor. Die Einträge in diesem Vektor sind auch hier std::tr1::array s, wobei die Länge der Arrays mit dem Index im Vektor zunimmt, daher sieht der Vektor wie folgt aus:

fusion::vector< std::tr1::array< double , 1 > , std::tr1::array< double , 2 > , std::tr1::array< double , 3 > , ... >

So vereinfacht sich dann beispielsweise die Koeffizienten-Definition des klassischen Runge-Kutta-Verfahrens zu:

using namespace std::tr1;
const array< double , 1 > coef_a1_rk4 = {{ 0.5 }};
const array< double , 2 > coef_a2_rk4 = {{ 0.0 , 0.5 }};
const array< double , 3 > coef_a3_rk4 = {{ 0.0 , 0.0 , 1.0 }};
const fusion::vector< array< double , 1 > , array< double , 2 > , array< double , 3 > > coef_a_rk4( coef_a1_rk4 , coef_a2_rk4 , coef_a3_rk4 );
const array< double , 4 > coef_b_rk4 = {{ 1.0 / 6.0 , 1.0 / 3.0 , 1.0 / 3.0 , 1.0 / 6.0 }};
const array< double , 4 > coef_c_rk4 = {{ 0.0 , 0.5 , 0.5 , 1.0 }};

Auch hier wird der eigentliche Algorithmus in mehrere Stages unterteilt. Diese sind exakt äquivalent zur MPL-Implementierung, außer dass die Koeffizienten nicht als Typen gespeichert werden, sondern als Array bzw. in Fusion-Vektoren. Der entscheidende Punkt ist dabei, dass diese Arrays/Vektoren als const deklariert werden und die Längen fix sind. Dadurch kann der Compiler (hoffentlich) Optimierungen anwenden, so als ob der Algorithmus direkt ausgeschrieben wird.

Die Rest der Implementierung ist dann:

template< class T , class Constant >
struct stage_fusion_wrapper
{
    typedef typename fusion::vector< size_t , T , std::tr1::array< T , Constant::value > > type;
};

template< class StateType , size_t NumOfStages >
class explicit_runge_kutta_fusion
{
    ...

    typedef typename fusion::result_of::as_vector
    <
        typename mpl::copy
        <
            mpl::range_c< int , 1 , num_of_stages + 1 > ,
            mpl::inserter
            <
                mpl::vector0<> ,
                mpl::push_back< mpl::_1 , stage_fusion_wrapper< double , mpl::_2 > >
            >
        >::type
    >::type stage_vector_base;

    struct stage_vector : public stage_vector_base
    {
        struct do_insertion
        {
            stage_vector_base &m_base;
            const coef_a_type &m_a;
            const coef_c_type &m_c;

            do_insertion( stage_vector_base &base , const coef_a_type &a , const coef_c_type &c )
            : m_base( base ) , m_a( a ) , m_c( c ) { }

            template< class Index >
            void operator()( Index ) const
            {
                fusion::at_c< 0 >( fusion::at< Index >( m_base ) ) = Index::value;
                fusion::at_c< 1 >( fusion::at< Index >( m_base ) ) = m_c[ Index::value ];
                fusion::at_c< 2 >( fusion::at< Index >( m_base ) ) = fusion::at< Index >( m_a );
            }
        };

        stage_vector( const coef_a_type &a , const coef_b_type &b , const coef_c_type &c )
        {
            typedef mpl::range_c< size_t , 0 , num_of_stages - 1 > indices;
            mpl::for_each< indices >( do_insertion( *this , a , c ) );
            fusion::at_c< 0 >( fusion::at_c< num_of_stages - 1 >( *this ) ) = num_of_stages - 1 ;
            fusion::at_c< 1 >( fusion::at_c< num_of_stages - 1 >( *this ) ) = c[ num_of_stages - 1 ];
            fusion::at_c< 2 >( fusion::at_c< num_of_stages - 1 >( *this ) ) = b;
        }
    };

    template< class System >
    struct calculate_stage
    {
        System &system;
        state_type &x , &x_tmp;
        state_type *k_vector;
        const double t;
        const double dt;

        calculate_stage( System &_system , state_type &_x , state_type &_x_tmp , state_type *_k_vector , const double _t , const double _dt )
        : system( _system ) , x( _x ) , x_tmp( _x_tmp ) , k_vector( _k_vector ) , t( _t ) , dt( _dt ) { }

        template< typename T , size_t stage_number >
        inline void operator()( fusion::vector< size_t , T , std::tr1::array< T , stage_number > > const &stage ) const
        {
            double c = fusion::at_c< 1 >( stage );
            if( stage_number == 1 )
                system( x , k_vector[stage_number-1] , t + c * dt );
            else
                system( x_tmp , k_vector[stage_number-1] , t + c * dt );

            if( stage_number == num_of_stages )
                fusion_algebra<stage_number>::foreach( x , x , fusion::at_c< 2 >( stage ) , k_vector , dt);
            else
                fusion_algebra<stage_number>::foreach( x_tmp , x , fusion::at_c< 2 >( stage ) , k_vector , dt);
        }
    };

public:

    explicit_runge_kutta_fusion( const coef_a_type &a , const coef_b_type &b , const coef_c_type &c )
    : m_a( a ) , m_b( b ) , m_c( c ) , m_stages( a , b , c ) , m_x_tmp() , m_k_vector() { }

    template< class System >
    inline void do_step( System &system , state_type &x , double t , const double dt )
    {
        fusion::for_each( m_stages , calculate_stage< System >( system , x , m_x_tmp , m_k_vector , t , dt ) );
    }

private:

    const coef_a_type m_a;
    const coef_b_type m_b;
    const coef_c_type m_c;
    const stage_vector m_stages;
    state_type m_x_tmp;
    state_type m_k_vector[num_of_stages];
};

Die einzelnen Stages werden hier auch wieder in einer eigenen Struktur gespeichert: stage_vector . Und calculate_stage übernimmt die Berechnung einer einzelnen Stage.

Performance

Die Performance der beiden Algorithmen ist etwa äquivalent. Ausserdem stehen beide Algorithmen den ausgeschriebenen Varianten in nichts hinterher, in einigen Fällen war der Fusion-Algorithmus sogar schneller als der Ausgeschriebene. Genaue Werte sind für verschiedene Compiler unterschiedlich, es gibt sogar Unterschiede zwischen gcc 4.4 und gcc 4.5.

Diskussion

Die Vorteile der beiden Ansätze liegen sofort auf der Hand. Man hat nur genau einen Algorithmus zu warten. Damit können praktisch alle klassischen expliziten Runge-Kutta Solver (wie Euler, RK4, Midpoint, Dormand-Prince 5-4, Dormand Prince 8-5-3, Fehlberg, ...) geschrieben werden. Ausserdem stehen diese Algorithmen den ausgeschriebenen Algorithmen in Performance in nichts nach. Man kann diese Konzept natürlich auch einfach auf die Stepper mit Fehlerberechnung oder auf Stepper mit kontinuierlichen Output erweitern.

Die Nachteile bei dieser Art von Programmierung sind gering. Der Code ist schwerer zu lesen als klassischer C/C++ Code und es wird auch relativ viel Boilerplate erzeugt. Allerdings wird die dadurch längere Entwicklungszeit durch die Allgemeinheit der Verfahren mehr als weggemacht.

Selbstverständlich erheben die beiden Version des Templatemetaprogrammierungs-Runge-Kutta-Verfahrens keinen Anspruch auf Allgemeingültigkeit. Sie wurde im Rahmen vom Boost.Odeint entwickelt und werden dort zum Einsatz kommen. Die Erweiterungsmöglichkeiten dazu sind vielzählig:

Abstrahierung des Zustandstyps. In den obigen beiden Beispielen wird der Zustand der Differentialgleichung durch std::tr1::array< double , N > beschrieben. Das kann selbstverständlich verallgemeinert werden.
Benutzungen von beliebigen Floating-Point Typen, wie float , std::complex< double > , GMP, ...
Erweiterung auf Solver mit Fehlerabschätzung und Dense-Output für Schrittweitensteuerung.
Verallgemeinerung der Algebra.
Höhere Ordnungen der Algebra. Im Moment werden nur Stepper bis Ordnung 4 unterstützt. Das ist trivial erweiterbar auf höhere Ordnungen. Auch können diese Erweiterungen durch Präprozessormakros automatisch erzeugt werden.

Über Kritik, Anregungen, Hinweise sowie Verbesserungsvorschläge würde ich mich sehr freuen.

Literatur

David Abrahams, Aleksey Gurtovoy, C++ Template Metaprogramming: Concepts, Tools, and Techniques from Boost and Beyond (Addison-Wesley, 2004).

Andrei Alexandrescu, Modern C++ Design, Generic Programming and Design Patterns Applied (Addison-Wesley, 2001)

Press William H et al., Numerical Recipes 3rd Edition: The Art of Scientific Computing, 3rd ed. (Cambridge University Press, 2007).

Ernst Hairer, Syvert P. Nørsett, and Gerhard Wanner, Solving Ordinary Differential Equations I: Nonstiff Problems, 2nd ed. (Springer, Berlin, 2009).

Milton Abramowitz, Irene A. Stegun, Handbook of Mathematical Functions (New York: Dover Publications, 1972). (online verfügbar)

GPC

So, eine Rechtschreibprüfung hab ich gemacht. War nicht viel zu tun, achte nur auf deine Apostrophen

Folgendermaßen geht's weiter: Wenn du keine Änderungen mehr am Artikel vornehmen willst, dann nimmst du die von mir gepostete Version und postest die hier in der Redaktion in einen neuen Thread mit dem [E] Tag. Sobald du dein Autorenprofil hast, verschiebe ich beides in die öffentlichen Foren.

headmyshoulder

GPC schrieb:

So, eine Rechtschreibprüfung hab ich gemacht. War nicht viel zu tun, achte nur auf deine Apostrophen

Folgendermaßen geht's weiter: Wenn du keine Änderungen mehr am Artikel vornehmen willst, dann nimmst du die von mir gepostete Version und postest die hier in der Redaktion in einen neuen Thread mit dem [E] Tag. Sobald du dein Autorenprofil hast, verschiebe ich beides in die öffentlichen Foren.

Ok, ich hab den Artikel so übernommen, und werd gleich den Artikel mit [E] posten. Ein Mini-Abschnitt mit dem Link zu den Sourcen ist noch dazugekommen. Das Autorenprofil hab ich auch in einem neuen Thread gepostet. Vielen Dank für deine nette Unterstützung bei der Arbeit an diesem Artikel.

Michael E.

Einführung

In diesem Artikel werde ich einige Techniken im Zusammenhang mit Templatemetaprogrammierung und numerischen Verfahren vorstellen. Insbesondere wird es dabei um numerische Verfahren zum Lösen von Differentialgleichungen (DGL) gehen. Die hier vorgestellten Methoden lassen sich aber problemlos auf eine Vielzahl von anderen numerischen Problemen anwenden.

Templatemetaprogrammierung ist eine moderne Programmiertechnik, bei welcher der C++ Compiler elementare Berechnungen durchführt oder komplexen C++-Code und Algorithmen erzeugt. Ein entscheidender Aspekt dabei ist, dass Templatemetaprogramme während des Kompiliervorganges ausgewertet bzw. ausgeführt werden. Das klassische Hello-World-Programm ist hier die Berechnung der Fakultät, welche direkt vom Compiler berechnet wird. Zur Laufzeit des fertigen Programmes wird diese Berechnung keine Zeit in Anspruch nehmen.

Templatemetaprogrammierung kann auch für komplexe Probleme und Algorithmen verwendet werden. Ein Beispiel ist die Bibliothek Boost.MetaStateMachine, mit welcher Zustandsmachinen anhand einer (statischen) Übergangstabelle während des Kompilierens erzeugt werden können. Ein weiteres Beispiel ist Boost.Accumulators, eine Bibliothek für statistische Auswertungen. Dort werden die genauen Berechnung und statistischen Routinen auch während der Kompilierphase erzeugt. Außerdem werden Techniken der Templatemetaprogrammierung auch oft im Zusammenhang mit Typkonvertierungen und Typüberprüfungen benutzt.

Es gibt eine Vielzahl von Bibliotheken, welche mehr oder weniger mit dem Thema Templatemetaprogrammierung zu tun haben. Die meisten davon sind in den Boost-Bibliotheken zu finden. Bevor man anfängt Metaprogramme selbst zu schreiben und zu implementieren, empfehle ich aber zuerst dort nachzusehen, ob nicht bereits eine adäquate Lösung existiert. Selbst durch das konsequente Anwenden und Benutzen dieser Bibliotheken kann man sehr viel lernen. Als wichtigste Bibliothek ist dabei sicherlich Boost.MPL mit einer Vielzahl von Metafunktionen, -algorithmen und Compile-Time-Sequenzen zu nennen~Satz war nicht vollständig~. Eine weitere wichtige Bibliothek ist Boost.Fusion, in welcher die Runtime-Aspekte mehr betont werden. Zusätzlich müssen in der Templatemetaprogrammierung sehr oft Typen klassifiziert werden. Dafür stehen etliche Funktion in Boost.TypeTraits zur Verfügung. Beispielsweise werden Metafunktionen wie is_pointer , is_class , aber auch Konversionen wie~fehlendes Wort~ add_const und remove_const zur Verfügung gestellt. Ein Teil dieser Funktionalität wird auch Einzug in den kommenden C++-Standard~Es ist ein Sprach- und kein Compiler-Standard~ finden.

Zur Entwicklung von Meta-Algorithmen werden wir hier ausschliesslich Boost.MPL verwenden. Diese Bibliothek ist in ihrem Umfang und ihren~fehlendes Wort~ Anwendungsmöglichkeiten am weitesten entwickelt. Es gibt schlicht und einfach keine wirklichen Alternativen dazu.

-----------------------------------------------

Wenn ich auch noch über den Rest gucken soll, müsst ihr noch ein paar Tage warten.

headmyshoulder

Michael E. schrieb:

Wenn ich auch noch über den Rest gucken soll, müsst ihr noch ein paar Tage warten.

Das wäre toll. Auf eine paar Tage mehr oder weniger kommt es nicht an. Danke

GPC

Michael E. schrieb:

Wenn ich auch noch über den Rest gucken soll, müsst ihr noch ein paar Tage warten.

Wäre klasse, aber dann nimm meine Version als Grundlage.

headmyshoulder

Michael E. schrieb:

Einführung

Wenn ich auch noch über den Rest gucken soll, müsst ihr noch ein paar Tage warten.

Wie siehts aus, wie lange wird es in etwas noch dauern?

Michael E.

Einführung

In diesem Artikel werde ich einige Techniken im Zusammenhang mit Templatemetaprogrammierung und numerischen Verfahren vorstellen. Insbesondere wird es dabei um numerische Verfahren zum Lösen von Differentialgleichungen (DGL) gehen. Die hier vorgestellten Methoden lassen sich aber problemlos auf eine Vielzahl von anderen numerischen Problemen anwenden.

Templatemetaprogrammierung ist eine moderne Programmiertechnik, bei welcher der C++ Compiler elementare Berechnungen durchführt oder komplexen C++-Code und Algorithmen erzeugt. Ein entscheidender Aspekt dabei ist, dass Templatemetaprogramme während des Kompiliervorganges ausgewertet bzw. ausgeführt werden. Das klassische Hello-World-Programm ist hier die Berechnung der Fakultät, welche direkt vom Compiler berechnet wird. Zur Laufzeit des fertigen Programmes wird diese Berechnung keine Zeit in Anspruch nehmen.

Templatemetaprogrammierung kann auch für komplexe Probleme und Algorithmen verwendet werden. Ein Beispiel ist die Bibliothek Boost.MetaStateMachine, mit welcher Zustandsmachinen anhand einer (statischen) Übergangstabelle während des Kompilierens erzeugt werden können. Ein weiteres Beispiel ist Boost.Accumulators, eine Bibliothek für statistische Auswertungen. Dort werden die genauen Berechnung und statistischen Routinen auch während der Kompilierphase erzeugt. Außerdem werden Techniken der Templatemetaprogrammierung auch oft im Zusammenhang mit Typkonvertierungen und Typüberprüfungen benutzt.

Es gibt eine Vielzahl von Bibliotheken, welche mehr oder weniger mit dem Thema Templatemetaprogrammierung zu tun haben. Die meisten davon sind in den Boost-Bibliotheken zu finden. Bevor man anfängt Metaprogramme selbst zu schreiben und zu implementieren, empfehle ich aber zuerst dort nachzusehen, ob nicht bereits eine adäquate Lösung existiert. Selbst durch das konsequente Anwenden und Benutzen dieser Bibliotheken kann man sehr viel lernen. Als wichtigste Bibliothek ist dabei sicherlich Boost.MPL mit einer Vielzahl von Metafunktionen, -algorithmen und Compile-Time-Sequenzen zu nennen~Satz war nicht vollständig~. Eine weitere wichtige Bibliothek ist Boost.Fusion, in welcher die Runtime-Aspekte mehr betont werden. Zusätzlich müssen in der Templatemetaprogrammierung sehr oft Typen klassifiziert werden. Dafür stehen etliche Funktion in Boost.TypeTraits zur Verfügung. Beispielsweise werden Metafunktionen wie is_pointer , is_class , aber auch Konversionen wie~fehlendes Wort~ add_const und remove_const zur Verfügung gestellt. Ein Teil dieser Funktionalität wird auch Einzug in den kommenden C++-Standard~Es ist ein Sprach- und kein Compiler-Standard~ finden.

Zur Entwicklung von Meta-Algorithmen werden wir hier ausschliesslich Boost.MPL verwenden. Diese Bibliothek ist in ihrem Umfang und ihren~fehlendes Wort~ Anwendungsmöglichkeiten am weitesten entwickelt. Es gibt schlicht und einfach keine wirklichen Alternativen dazu.

Wer sich ernsthaft mit dem Thema Templatemetaprogrammierung auseinandersetzen will, dem würde ich einige Pflichtlektüren empfehlen. "Modern C++ Design" von Andrei Alexandrescu ist vielleicht ein guter, wenn auch nicht einfacher Einstieg. Dort werden so gut wie alle der grundlegenden Techniken erklärt, aber man sollte schon recht sicher im Umgang mit C++ sein. Das Buch "C++ Template Metaprogramming" von David Abrahams und Aleksey Gurtovoy erklärt die Thematik recht anschaulich am Beispiel von Boost.MPL, ohne allerdings detailliert die Implementierungen der einzelnen Algorithemen zu erläutern. Und zu guter Letzt sei noch das Buch "C++ Templates" von David Vandevoorde und Nicolai Josuttis erwähnt, welches sich ausführlichst mit Templates beschäftigt und damit die (theoretischen) Grundlagen für vielen Methoden liefert.

Der Artikel ist wie folgt gegliedert. Nach dieser Einführung werde ich im zweiten Teil zeigen, wie numerische Konstanten für Templatemetaprogrammierung definiert werden. Wir werden diese Definition später, wenn es um die Konstruktion der numerische Algorithmen geht, wieder benötigen. Danach werde ich anhand eines Beispieles kurz erläutern, welche Probleme generell beim Erstellen von Templatemetaprogrammen auftreten können, wie man am besten Fehler in diesen Programmen findet und Templatemetaprogramme debuggt. Danach kommt der Hauptteil, in welchem die klassischen Runge-Kutta-Verfahren zum Lösen von Differentialgleichungen mit Hilfe von Templatemetaprogrammen konstruiert werden. Der Artikel endet mit einer Diskussion der Ergebnisse und Methoden.

Numerische Konstanten

Numerische Konstanten können mit den eingebauten integralen Typen der MPL dargestellt werden:

typedef mpl::int_< 0 > zero;
typedef mpl::char_< 1 > one;
typedef mpl::long_< 2 > two;
typedef mpl::long_< -2 > minus_two;

Der numerische Wert wird gewöhnlich in einer statischen Membervariable ::value abgelegt und kann wie folgt ausgegeben werden.

cout << zero::value << endl;
cout << one::value << endl;
cout << two::value << endl;
cout << minus_two::value << endl;

Ausserdem sind diese Typen nullary Metafunktionen, welche ihren Typ zurückgeben. Daher sind die beiden Typen mpl::int_< 1 > und mpl::int_< 1 >::type gleichwertig.

Man beachte: Da Zahlen als Typen definiert sind, steht ihr Wert während des Kompiliervorgangs fest. Ausdrücke, welche diese Werte benutzen, können dann vom Kompiler optimiert werden. Allerdings können nur die ganzzahligen Typen (***char, short, int, long, ...***) mit Hilfe der numerischen Konstanten der MPL ausgedrückt werden. Benötigen wir Fließkommazahlen zur Compilezeit, so können wir diese durch rationale Zahlen darstellen:

template< class T , T Numerator , T Denominator >
struct rational_number
{
    typedef T value_type;
    typedef rational_number type;
    const static value_type numerator = Numerator;
    const static value_type denominator = Denominator;
};

template< class Num >
struct conv
{
    static double get_value( void )
    {
        return static_cast< double >( Num::value );
    }
};

template< class T , T Numerator , T Denominator >
struct conv< rational_number< T , Numerator , Denominator > >
{
    typedef rational_number< T , Numerator , Denominator > num_type;
    static double get_value( void )
    {
        return static_cast< double >( num_type::numerator ) / static_cast< double >( num_type::denominator );
    }
};

Die Klasse conv übernimmt hier nur die Konvertierung zu double s. Man sieht aber hier auch einen bedeutenden Unterschied zu integralen Typen. Der Wert einer rationalen Zahl kann nicht so ohne Weiteres in einer statischen Member-Variable abgelegt werden, da dies nur für integrale Typen erlaubt ist, es sein denn, man definiert diese Variable extern, was aber nicht besonders ansprechend ist. Außerdem ist die Variable dann auch wirklich als Symbol vorhanden. Daher muss hier der Umweg über eine statische Klasse gegangen werden. Diese Klasse ist spezialisiert für rational_number und man kann mit ihr sowohl integrale als auch Fließkommazahlen einheitlich benutzen, beispielsweise so:

typedef rational_number< long , 10 , 15 > two_third;
cout << conv< two_third >::get_value() << endl;
cout << conv< two >::get_value() << endl;

Für unsere Anwendung ist die Definition von numerischen Konstanten völlig ausreichend. Zusätzlich könnte man jetzt noch Operation *(+, -, , /) definieren. Es gibt sogar einige Bibliotheken, welche die Benutzung von numerischen Konstanten soweit auf die Spitze treiben, dass beispielsweise trigonometrische Funktionen wie der Sinus oder der Cosinus vom Compiler berechnet werden können, siehe hier.

Beispielalgorithmen und Debugging

Debugging von Templatemetaalgorithmen ist meistens problematisch. Eine richtige Lösung hab ich nicht, allerdings kann man mit Hilfe von typeid einige Informationen über mögliche Probleme erhalten. Nehmen wir ein Beispiel: Wir wollen einen beliebig langen [c]fusion[/c]-Vektor konstruieren, dessen Elemente tr1::array s sind. Die Länge der Arrays soll dabei ihrem Index in dem Vektor entsprechen. Beispielsweise soll ein Vektor der Länge 3 so aussehen

typedef fusion::vector
<
    tr1::array< double , 1 > ,
    tr1::array< double , 2 > ,
    tr1::array< double , 3 >
> coef_type;

Ein erster Versuch solch einen Vektor zu konstruieren wäre

typedef mpl::copy
<
    mpl::range_c< size_t , 1 , N + 1 > ,
    mpl::inserter
    <
        mpl::vector<> , 
        mpl::push_back
        <
            mpl::_1 ,
            tr1::array< double , mpl::_2::value >
        >
    >
>::type mpl_coef_type;
typedef fusion::result_of::as_vector< mpl_coef_type >::type coef_type;

Dies würde einem klassischen C++ Kopieralgorithmus mit einem speziellen back_insert_iterator entsprechen, ganz analog zu

// Achtung, das Beispiel wird nicht kompilieren
vector< size_t > indices;
for( size_t i=1 ; i<=N ; ++i ) indices.push_back( i );
vector< arrays > vec;
copy( indices.begin() , indices.begin() , back_insert_iterator< vector< arrays > >( vec ) );

Aber gehen wir Schritt für Schritt durch den Meta-Algorithmus:

mpl::copy< Seq , Inserter > ist der zentrale Kopiervorgang. Er nimmt eine MPL-Sequenz Seq entgegen und gibt jedes Element dieser Sequenz an den Inserter weiter.
mpl::range_c< size_t , 1 , N + 1 > konstruiert einen [c]mpl::vector[/c], dessen Elemente size_t s von 1 bis N sind, also wie mpl::vector< 1 , 2 , 3 , ... , N-1 , N > . Dieser Vektor wird benötigt, da er erstens die richtige Anzahl von Elementen für den finalen Vektor besitzt und zweitens jedes Element die Länge des aktuellen Arrays ist.
mpl::inserter< Seq, Op > ist ein allgemeiner Inserter, ähnlich zu den Insert-Iteratoren der STL. Als erstes Template-Argument wird hier eine initiale Sequenz erwartet, in unserem Falle ein leerer Vektor mpl::vector<> . Das zweite Argument definiert die Einfüge-Operation.
mpl::push_back< Seq , El > macht nichts anderes als El an die Sequenz Seq anzuhängen. In unserem Fall wird also an den aktuellen Vektor ein neues Element vom Typ std::tr1::array angehängt. Hier sieht man auch sehr schön die Verwendung von Platzhaltern. mpl::_1 und mpl::_2 bezeichnen das erste bzw. zweite Argument, welches an push_back übergeben wird.
fusion::result_of::as_vector konvertiert den Mpl-Vektor in einen Fusion-Vektor.

Ausgeschrieben soll wird der gesamte Algorithmus in etwa so aussehen~Ich verstehe den Satz nicht~:

mpl::push_back
<
    mpl::push_back
    < 
        mpl::push_back
        <
            mpl::vector<> ,
            tr1::array< double , 1 >
        >::type ,
        tr1::array< double , 2 >
    >::type ,
    tr1::array< double , 3 >
>::type

Der Code kompiliert und man könnte meinen, dass wir die Aufgabe geschafft haben. Aber dem ist nicht so. Wollen wir diesen jetzt verwenden, wird der Compiler eine schier unglaubliche Anzahl von Fehlern ausspucken. Diese nach ihrer Essenz zu durchforsten, ist mühsam und selbst wenn man die entscheidenden Stellen gefunden hat, kann es sein, dass man mit der Fehlermeldung nicht viel anfangen kann (mir geht das zumindest so). Um diese Fehlersuche abzukürzen, gibt es mehrere Wege:

Benutze STLFilt! Es wird die Anzahl der Fehlermeldungen drastisch reduzieren.
Benutze statische Assertions! Damit können bestimmte Fehler viel früher gefunden werden. Statische Assertions sind insbesondere für Library-Entwickler wichtig, weil damit Bedingungen an die Typen geprüft werden.
Benutze einen Typdebugger! Er wird Information über den Typ liefern, welche eventuell hilfreich sind.

template< class T >
struct wrap {};

struct debug_type
{
    size_t m_indent;

    debug_type( size_t indent = 0 )
    : m_indent( indent ) { }

    template< class T >
    void operator()( wrap<T> ) const
    {
        typedef typename boost::mpl::is_sequence< T >::type is_seq;
        debug_impl( wrap< T >() , is_seq() );
    }

    template< class Seq >
    void debug_impl( wrap< Seq > , boost::mpl::true_ ) const
    {
        std::cout << indent() << "Sequence : " << typeid(Seq).name() << std::endl;
        boost::mpl::for_each<
        typename boost::mpl::transform< Seq , wrap< boost::mpl::_1 > >::type >( debug_type( m_indent + 1 ) );
    }

    template< class T >
    void debug_impl( wrap< T > , boost::mpl::false_ ) const
    {
        std::cout << indent() << typeid(T).name() << std::endl;
    }

    std::string indent( void ) const
    {
        std::string in;
        for( size_t i=0 ; i<m_indent ; ++i ) in += "  ";
        return in;
    }
};

template< class T >
void debug( size_t indent = 0 )
{
    debug_type d( indent );
    wrap< T > t;
    d( t );
}

Der wrap -Trick kommt aus dem Buch von David Abrahams und sorgt dafür, dass der Typ-Printer auch mit Referenzen funktioniert (Referenzen sind nicht default-konstruierbar und können nicht einfach von mpl::for_each erzeugt werden).

In unserem Beispiel des Vektors von Arrays führt nun die Benutzung von

debug< mpl_coef_type >();
debug< coef_type >();
dout << endl;

zu folgender Ausgabe:

Sequence : N5boost3mpl6v_itemINSt3tr15arrayIdLm2EEENS1_IS4_NS1_IS4_NS0_6vectorIN4mpl_2naES7_S7_S7_S7_S7_S7_S7_S7_S7_S7_S7_S7_S7_S7_S7_S7_S7_S7_S7_EELi0EEELi0EEELi0EEE
  NSt3tr15arrayIdLm2EEE
  NSt3tr15arrayIdLm2EEE
  NSt3tr15arrayIdLm2EEE
Sequence : N5boost6fusion7vector3INSt3tr15arrayIdLm2EEES4_S4_EE
  NSt3tr15arrayIdLm2EEE
  NSt3tr15arrayIdLm2EEE
  NSt3tr15arrayIdLm2EEE

Und voilà, man sieht, dass mpl::_2::value nicht richtig aufgelöst/ersetzt wurde, denn die Arrays haben die Länge 2. Eine Lösung ist schnell gefunden, wir packen die Definition des Arrays einfach in eine Meta-Funktion:

template< class T , class Constant >
struct array_wrapper
{
    typedef typename tr1::array< T , Constant::value > type;
};

typedef mpl::copy
<
    mpl::range_c< size_t , 1 , N + 1 > ,
    mpl::inserter
    <
        mpl::vector<> ,
        mpl::push_back
        <
            mpl::_1 ,
            array_wrapper< double , mpl::_2 >
        >
    >
>::type mpl_coef_type;
typedef fusion::result_of::as_vector< mpl_coef_type >::type coef_type;

und unser kleiner, aber feiner Algorithmus liefert das gewünschte Ergebnis.

Konstruktion von numerischen Algorithmen

Kommen wir nun zum eigentlich Thema dieses Artikels, der Konstruktion von numerischen Algorithmen mit Hilfe von Template-Metaprogramming. Wir verdeutlichen unseren Ansatz anhand der klassischen expliziten Runge-Kutta-Solvern zum Lösen von gewöhnlichen Differentialgleichungen. Ähnliche Ansätze können auch für verwandte numerische Probleme wie

numerische Integration / Lagrange-Integration,
Konstruktion von Spline-Basis-Funktionen oder auch
finite Differenzen zum Lösen von partiellen Differentialgleichungen

genutzt werden.

Ein gewöhnliche Differentialgleichung beschreibt die (zeitliche) Entwicklung einer Größe x anhand ihrer Ableitung(en):

dx / dt = f(x,t)

x kann dabei vektoriell sein, daher aus mehreren Werten bestehen. Existiert zu einem Zeitpunkt t=t0 eine Anfangsbedingung x(t0)=x0 und ist die zeitliche Entwicklung von x gesucht, so spricht man von einem Anfangswertproblem. Numerisch werden Anfangswertprobleme normalerweise iterativ gelöst. Man bestimmt also eine aufeinanderfolgende Sequenz x(t0),x(t0+dt),x(t0+2dt),... , welche die Lösung numerisch repäsentiert. dt ist dabei die Schrittweite.

Eine einfache Klasse von solchen Algorithmen sind sogenannte explizite Runge-Kutta-Verfahren. Sie haben den großen Vorteil, dass keine Gleichungssysteme gelöst werden müssen. Schematisch werden Runge-Kutta-Solver durch Butcher-Tableaus dargestellt:

c[t]1[/t] |
c[t]2[/t] | a[t]2,1[/t]
c[t]3[/t] | a[t]3,1[/t] a[t]3,2[/t]
...
c[t]s[/t] | a[t]s,1[/t] a[t]s,2[/t] ... a[t]s,s-1[/t]
-----------------------------------
   | b[t]1[/t] b[t]2[/t] ... b[t]s-1[/t] b[t]s[/t]

s bezeichnet dabei die Ordnung des Verfahrens, was ein wichtiger Parameter ist, da er die Genauigkeit der Lösung charakterisiert.

Aus solch einer Tabelle kann man nun das Verfahren zum Lösen der DGL wie folgt konstruieren:

k[t]1[/t] = f( x(t) , t + c[t]1[/t] * dt )
k[t]2[/t] = f( x(t) + dt * a[t]2,1[/t] * k[t]1[/t] , t + c[t]2[/t] * dt )
k[t]3[/t] = f( x(t) + dt * ( a[t]3,1[/t] * k[t]1[/t] + a[t]3,2[/t] * k[t]2[/t] ), t + c[t]3[/t] * dt )
...
k[t]s[/t] = f( x(t) + dt * ( a[t]s,1[/t] * k[t]1[/t] + ... + a[t]s,s-1[/t] * k[t]s-1[/t] ) , t + c[t]s[/t] * dt )
x(t+dt) = x(t) + dt * ( b[t]1[/t] k[t]1[/t] + ... + b[t]s[/t] k[t]s[/t] )

Ein einfaches Beispiel ist das Euler-Verfahren:

0 |
--|---
  | 1

Oder ausgeschrieben:

k[t]1[/t] = f( x(t) , t )
x(t+dt) = x(t) + dt * k[t]1[/t]

Ein etwas komplizierteres Beispiel ist das klassische Runge-Kutta-Verfahren vierter Ordnung (im Folgenden als RK4 bezeichnet):

0   |
1/2 | 1/2
1/2 | 0   1/2
1   | 0   0   1
----|----------------
    | 1/3 1/6 1/6 1/3

In Pseudocode sieht diese Methode so aus:

k[t]1[/t] = f( x(t) , t )
k[t]2[/t] = f( x(t) + dt * 1/2 * k[t]1[/t] , t + 1/2 * dt )
k[t]3[/t] = f( x(t) + dt * 1/2 * k[t]2[/t] , t + 1/2 * dt )
k[t]4[/t] = f( x(t) + dt * 1 * k[t]3[/t] , t + 1 * dt )
x(t+dt) = x(t) + dt * ( 1/6 k[t]1[/t] + 1/3 k[t]2[/t] + 1/3 k[t]3[/t] + 1/6 k[t]4[/t] )

Wir werden nun Schritt für Schritt einen Meta-Algorithmus konstruieren, welcher solche Algorithmen während des Kompiliervorganges erzeugt. Als erstes werden die Koeffizienten a[t]i,j[/t] , b[t]i[/t] und c[t]i[/t] definiert. Danach wird das Hauptgerüst des Algorithmus entworfen und im letzten Teil die einzelnen Schritte ausgeführt. Zum Schluss wird eine alternative Implementation vorgestellt, welche sich durch eine einfache Konstrukion der Koeffizienten auszeichnet.

Die Koeffizienten eines Verfahrens werden in einem mpl::vector abgelegt. Für RK4 sieht die Definition wie folgt aus:

// allgemeine Konstanten
typedef mpl::int_< 0 > null;
typedef mpl::int_< 1 > one;
typedef rational_number< long , 1 , 2 > one_half;
typedef rational_number< long , 1 , 3 > one_third;
typedef rational_number< long , 1 , 6 > one_sixth;

typedef mpl::vector
<
    mpl::vector< one_half > ,
    mpl::vector< null , one_half > ,
    mpl::vector< null , null , one >
> rk4_a;
typedef mpl::vector< null , one_half , one_half , one > rk4_c;
typedef mpl::vector< one_sixth , one_third , one_third , one_sixth > rk4_b;

Der Solver selbst ist eine templatisierte Klasse, welche durch die Koeffizienten parametrisiert ist:

template< class State , class ButcherA , class ButcherB , class ButcherC >
class explicit_runge_kutta
{
public:
    typedef ButcherA butcher_a;
    typedef ButcherB butcher_b;
    typedef ButcherC butcher_c;

    static const size_t num_of_sages = mpl::size< ButcherC >::value;

    typedef double time_type;
    typedef State state_type;
...
}

time_type und state_type geben die Typen von t bzw. x an und für state_type wird ein std::tr1::array< double , N > erwartet.

Als Nächstes überprüfen wir, ob butcher_a , butcher_b und butcher_c den Anforderungen an die Koeffizienten genügen. Insbesondere prüfen wir, ob diese Typlisten die entsprechenden Anzahlen von Elementen enthalten.

template< class State , class ButcherA , class ButcherB , class ButcherC >
class explicit_runge_kutta
{
    ...
    BOOST_STATIC_ASSERT( ( num_of_stages > 0 ) );
    BOOST_STATIC_ASSERT( ( num_of_stages == size_t( mpl::size< butcher_b >::value ) ) );
    BOOST_STATIC_ASSERT( ( ( num_of_stages - 1 ) == size_t( mpl::size< butcher_a >::value ) ) ) ;

    typedef typename mpl::equal
    <
        mpl::range_c< int , 1 , num_of_stages > ,
        typename mpl::transform< butcher_a , mpl::size< mpl::_1 > >::type ,
        mpl::equal_to< mpl::_1 , mpl::_2 >
    >::type equal_type;
    BOOST_STATIC_ASSERT(( equal_type::value == true ));
};

Jetzt kommen wir zum eigentlichen Kernstück unseres Metaalgorithmus. Die expliziten Runge-Kutta-Verfahren verfügen über s Stufen, wobei s hier wieder die Ordnung des Verfahrens ist. Die erste Stufe sieht wie folgt aus:

k[t]1[/t] = f( x(t) , t + c[t]1[/t] )
xtmp = x(t) + dt * a[t]2,1[/t] * k[t]1[/t]

In den nachfolgenden Stufen (außer der letzten) wird immer

k[t]i[/t] = f( xtmp , t + c[t]i[/t] )
xtmp = x(t) + dt * a[t]i+1,1[/t] * k[t]1[/t] + ... + a[t]i+1,i[/t] * k[t]i[/t]

berechnet. Die letzte Stufe unterscheidet sich von den anderen~Ich hoffe, ich habe den Sinn getroffen, das "auch" verstehe ich nicht~, hier wird das Ergebnis geschrieben:

k[t]s[/t] = f( xtmp , t + c[t]s[/t] )
x(t+dt) = x(t) + dt * b[t]1[/t] * k[t]1[/t] + ... + b[t]s[/t] * k[t]s[/t]

Diesen Algorithmus setzten wir nun in einen Meta-Algorithmus um:

template< class State , class ButcherA , class ButcherB , class ButcherC >
class explicit_runge_kutta_mpl
{
public:
    ...
    typedef mpl::range_c< size_t , 0 , num_of_stages > indices;
    typedef typename mpl::push_back< butcher_a , butcher_b >::type butcher_right;
    typedef mpl::zip_view< mpl::vector< indices , butcher_c , butcher_right > > butcher_tableau;

    template< class System >
    void do_step( System &system , state_type &x , double t , const double dt )
    {
        mpl::for_each< butcher_tableau >( calculate_stage< System >( system , x , m_x_tmp , m_k_vector , t , dt ) );
    }
private:

    state_type m_k_vector[num_of_stages];
    state_type m_x_tmp;
    ...
};

Schauen wir uns an, was hier passiert. butcher_tableau ist ein Mpl-Vektor mit folgenden Einträgen:

0 , c[t]1[/t] , a[t]2,1[/t] 
1 , c[t]2[/t] , a[t]3,1[/t] , a[t]3,2[/t]
2 , c[t]3[/t] , a[t]4,1[/t] , a[t]4,2[/t] , a[t]4,3[/t] 
...
s , c[t]s[/t] , b[t]1[/t]   , b[t]2[/t]   , b[t]3[/t]   , ... , b[t]s[/t]

Die i-te Zeile enthält also alle Informationen, welche wir zur Berechnung der i-ten Stufe benötigen. Über diesen Vektor wird nun in der Zeile

mpl::for_each< butcher_tableau >( calculate_stage< System >( system , x , m_x_tmp , m_k_vector , t , dt ) );

iteriert und in jeder Iteration calculate_state aufgerufen. Die Methode do_step geht nun einen Zeitschritt in der Differentialgleichung nach vorn, do_step berechnet x(t+dt) anhand von x(t) . Die Differentialgleichung selbst wird durch system angegeben. system kann hier ein Funktor, aber auch ein Funktionspointer sein. Es wird dabei immer erwartet, dass wir

state_type dxdt;
system( x , dxdt , t );

ausführen können.

Nun wenden wir uns der der eigentlichen Berechnung jeder einzelnen Stufe zu. Das passiert in calculate_stage , was eine Member-Struktur von explicit_runge_kutta_mpl ist. Diese Struktur ist gleichzeitig auch ein Funktor, daher existiert der Klammer-Operator.

template< class State , class ButcherA , class ButcherB , class ButcherC >
class explicit_runge_kutta_mpl
{
public:
    ...
    template< class System >
    struct calculate_stage
    {
        System &system;
        state_type &x , &x_tmp;
        state_type *k_vector;
        const double t;
        const double dt;

        calculate_stage( System &_system , state_type &_x , state_type &_x_tmp , state_type *_k_vector , const double _t , const double _dt )
        : system( _system ) , x( _x ) , x_tmp( _x_tmp ) , k_vector( _k_vector ) , t( _t ) , dt( _dt )
        {}

        template< class Stage >
        void operator()( Stage v )
        {
            typedef typename mpl::at< Stage , mpl::int_< 0 > >::type index_type;
            typedef typename mpl::at< Stage , mpl::int_< 1 > >::type time_value_type;
            typedef typename mpl::at< Stage , mpl::int_< 2 > >::type coef_type;

            const static size_t index = index_type::value;
            const static double time_value = conv< time_value_type >::get_value();

            BOOST_STATIC_ASSERT(( ( mpl::size< coef_type >::value - 1 ) == index ));

            if( index == 0 ) system( x , k_vector[index] , t + time_value * dt );
            else system( x_tmp , k_vector[index] , t + time_value * dt );

            if( index == ( num_of_stages - 1 ) )
                algebra< state_type , coef_type , mpl::int_< index > >::do_step( x , x , k_vector , dt );
            else
                algebra< state_type , coef_type , mpl::int_< index > >::do_step( x_tmp , x , k_vector , dt );
        }
    };
    ...
};

Das Klassen-Template algebra enthält nur eine einzelne statische Methode do_step und ist durch den Zustandstyp state_type , die Koeffizienten des Verfahrens coef_type und die Ordnung des Verfahrens order parametrisiert.

template< class state_type , class coef_tye , class order > struct algebra_mpl;

algebra muss für alle Ordnungen spezialisiert werden. Hier könnte bestimmt der Präprozessor helfen, aber das ist eine andere Geschichte. Deswegen schreiben wir alle Ordnungen aus. Zum Beispiel sehen die ersten beiden Ordnungen so aus:

template< class state_type , class coef_type >
struct algebra_mpl< state_type , coef_type , mpl::int_< 0 > >
{
    static void do_step( state_type &x_tmp , const state_type &x , const state_type *k_vector , const double dt )
    {
        const static double a1 = dt * conv< typename mpl::at< coef_type , mpl::int_< 0 > >::type >::get_value();

        for( size_t i=0 ; i<x.size() ; ++i )
        {
            x_tmp[i] = x[i] + a1 * k_vector[0][i];
        }
    }
};

template< class state_type , class coef_type >
struct algebra_mpl< state_type , coef_type , mpl::int_< 1 > >
{
    static void do_step( state_type &x_tmp , const state_type &x , const state_type *k_vector , const double dt )
    {
        const static double a1 = dt * conv< typename mpl::at< coef_type , mpl::int_< 0 > >::type >::get_value();
        const static double a2 = dt * conv< typename mpl::at< coef_type , mpl::int_< 1 > >::type >::get_value();

        for( size_t i=0 ; i<x.size() ; ++i )
        {
            x_tmp[i] = x[i] + a1 * k_vector[0][i] + a2 * k_vector[1][i];
        }
    }
};

Damit sind wir fertig. Das Letzte, was uns jetzt noch übrig bleibt, ist einen speziellen Solver zu initialisieren sowie zu instanzieren und dann können wir mit dem Lösen von Differentialgleichungen beginnen. Für das vorangegangene Euler-Beispiel sieht die Definition des Solvers wie folgt aus:

typedef mpl::int_< 0 > null;
typedef mpl::int_< 1 > one;

typedef mpl::vector< > euler_a;
typedef mpl::vector< one > euler_b;
typedef mpl::vector< null > euler_c;

template< class State > class stepper_euler
: public explicit_runge_kutta_mpl< State , euler_a , euler_b , euler_c > {};

Als erstes werden die Konstanten für das Verfahren definiert, in diesem Fall 0 und 1. Danach werden die Koeffizienten a[t]i,j[/t] , b[t]i[/t] und c[t]i[/t] definiert. Die Koeffizienten b und c enthalten jeweils nur einen Eintrag, die Koeffizienten a sind sogar leer. Danach wird der Euler-Solver per Vererbung definiert. Für den klassischen Runge-Kutta-Solver sieht alles sehr ähnlich aus, nur dass andere Koeffizienten vorkommen:

typedef rational_number< int , 1 , 2 > one_half;
typedef rational_number< int , 1 , 3 > one_third;
typedef rational_number< int , 1 , 6 > one_sixth;

typedef mpl::vector
<
    mpl::vector< one_half > ,
    mpl::vector< null , one_half > ,
    mpl::vector< null , null , one >
> rk4_a;
typedef mpl::vector< null , one_half , one_half , one > rk4_c;
typedef mpl::vector< one_sixth , one_third , one_third , one_sixth > rk4_b;

template< class State > class stepper_rk4
: public explicit_runge_kutta_mpl< State , rk4_a , rk4_b , rk4_c > { };

Alternative [kor]Implementierung[/kor]

Ein Nachteil der obigen Implementierung des Compilezeit-Runge-Kutta Algorithmus ist sicherlich die Art und Weise, wie die Koeffizienten angegeben werden müssen. Für einen Koeffizienten wie 1/6 müssen wir beispielsweise rational_number< int , 1 , 6 > schreiben. Eine einfachere Definition könnte beispielsweise konstante Arrays benutzen, welche im Konstruktor des Runge-Kutta-Verfahrens übergeben werden. Wir werden diesen Ansatz jetzt entwickeln, er erfordert allerdings ein kleines bisschen mehr Metaprogrammierungszauberei.

In der alternativen Implementierung werden die Typen der Koeffizienten folgendermaßen definiert:

template< class T , class Constant >
struct array_wrapper
{
    typedef typename std::tr1::array< T , Constant::value > type;
};

template< class StateType , size_t NumOfStages >
class explicit_runge_kutta_fusion
{
public:
    typedef StateType state_type;
    const static size_t num_of_stages = NumOfStages;

    typedef typename fusion::result_of::as_vector
    <
        typename mpl::copy
        <
            mpl::range_c< size_t , 1 , num_of_stages > ,
            mpl::inserter
            <
                mpl::vector0< > ,
                mpl::push_back< mpl::_1 , array_wrapper< double , mpl::_2 > >
            >
        >::type
    >::type coef_a_type;

    typedef std::tr1::array< double , num_of_stages > coef_b_type;
    typedef std::tr1::array< double , num_of_stages > coef_c_type;
    ...
}

Die Koeffizienten b und c werden in einfachen Arrays gespeichert. Für die Koeffizienten a benötigen wir allerdings etwas Komplizierteres: einen fusion::vector . Solche Vektoren sind ähnlich zu std::tr1::tuple , aber mit dem Unterschied, dass dazu eine riesige Menge an Meta-Algorithmen und Funktionen zur Verfügung stehen. Der Meta-Algorithmus hier konstruiert einen solchen Vektor (siehe auch das obige Beispiel). result_of::as_vector konvertiert einen MPL-Vektor in einen Fusion-Vektor. Die Einträge in diesem Vektor sind auch hier std::tr1::array s, wobei die Länge der Arrays mit dem Index im Vektor zunimmt, daher sieht der Vektor wie folgt aus:

fusion::vector< std::tr1::array< double , 1 > , std::tr1::array< double , 2 > , std::tr1::array< double , 3 > , ... >

So vereinfacht sich dann beispielsweise die Koeffizienten-Definition des klassischen Runge-Kutta-Verfahrens zu

using namespace std::tr1;
const array< double , 1 > coef_a1_rk4 = {{ 0.5 }};
const array< double , 2 > coef_a2_rk4 = {{ 0.0 , 0.5 }};
const array< double , 3 > coef_a3_rk4 = {{ 0.0 , 0.0 , 1.0 }};
const fusion::vector< array< double , 1 > , array< double , 2 > , array< double , 3 > > coef_a_rk4( coef_a1_rk4 , coef_a2_rk4 , coef_a3_rk4 );
const array< double , 4 > coef_b_rk4 = {{ 1.0 / 6.0 , 1.0 / 3.0 , 1.0 / 3.0 , 1.0 / 6.0 }};
const array< double , 4 > coef_c_rk4 = {{ 0.0 , 0.5 , 0.5 , 1.0 }};

Auch hier wird der eigentliche Algorithmus in mehrere Stages unterteilt. Diese sind exakt äquivalent zur MPL-Implementierung, außer dass die Koeffizienten nicht als Typen gespeichert werden, sondern als Array bzw. in Fusion-Vektoren. Der entscheidende Punkt ist dabei, dass diese Arrays/Vektoren als const deklariert werden und die Längen fix sind. Dadurch kann der Compiler (hoffentlich) Optimierungen anwenden, so als ob der Algorithmus direkt ausgeschrieben würde.

Der Rest der Implementierung ist dann:

template< class T , class Constant >
struct stage_fusion_wrapper
{
    typedef typename fusion::vector< size_t , T , std::tr1::array< T , Constant::value > > type;
};

template< class StateType , size_t NumOfStages >
class explicit_runge_kutta_fusion
{
    ...

    typedef typename fusion::result_of::as_vector
    <
        typename mpl::copy
        <
            mpl::range_c< int , 1 , num_of_stages + 1 > ,
            mpl::inserter
            <
                mpl::vector0<> ,
                mpl::push_back< mpl::_1 , stage_fusion_wrapper< double , mpl::_2 > >
            >
        >::type
    >::type stage_vector_base;

    struct stage_vector : public stage_vector_base
    {
        struct do_insertion
        {
            stage_vector_base &m_base;
            const coef_a_type &m_a;
            const coef_c_type &m_c;

            do_insertion( stage_vector_base &base , const coef_a_type &a , const coef_c_type &c )
            : m_base( base ) , m_a( a ) , m_c( c ) { }

            template< class Index >
            void operator()( Index ) const
            {
                fusion::at_c< 0 >( fusion::at< Index >( m_base ) ) = Index::value;
                fusion::at_c< 1 >( fusion::at< Index >( m_base ) ) = m_c[ Index::value ];
                fusion::at_c< 2 >( fusion::at< Index >( m_base ) ) = fusion::at< Index >( m_a );
            }
        };

        stage_vector( const coef_a_type &a , const coef_b_type &b , const coef_c_type &c )
        {
            typedef mpl::range_c< size_t , 0 , num_of_stages - 1 > indices;
            mpl::for_each< indices >( do_insertion( *this , a , c ) );
            fusion::at_c< 0 >( fusion::at_c< num_of_stages - 1 >( *this ) ) = num_of_stages - 1 ;
            fusion::at_c< 1 >( fusion::at_c< num_of_stages - 1 >( *this ) ) = c[ num_of_stages - 1 ];
            fusion::at_c< 2 >( fusion::at_c< num_of_stages - 1 >( *this ) ) = b;
        }
    };

    template< class System >
    struct calculate_stage
    {
        System &system;
        state_type &x , &x_tmp;
        state_type *k_vector;
        const double t;
        const double dt;

        calculate_stage( System &_system , state_type &_x , state_type &_x_tmp , state_type *_k_vector , const double _t , const double _dt )
        : system( _system ) , x( _x ) , x_tmp( _x_tmp ) , k_vector( _k_vector ) , t( _t ) , dt( _dt ) { }

        template< typename T , size_t stage_number >
        inline void operator()( fusion::vector< size_t , T , std::tr1::array< T , stage_number > > const &stage ) const
        {
            double c = fusion::at_c< 1 >( stage );
            if( stage_number == 1 )
                system( x , k_vector[stage_number-1] , t + c * dt );
            else
                system( x_tmp , k_vector[stage_number-1] , t + c * dt );

            if( stage_number == num_of_stages )
                fusion_algebra<stage_number>::foreach( x , x , fusion::at_c< 2 >( stage ) , k_vector , dt);
            else
                fusion_algebra<stage_number>::foreach( x_tmp , x , fusion::at_c< 2 >( stage ) , k_vector , dt);
        }
    };

public:

    explicit_runge_kutta_fusion( const coef_a_type &a , const coef_b_type &b , const coef_c_type &c )
    : m_a( a ) , m_b( b ) , m_c( c ) , m_stages( a , b , c ) , m_x_tmp() , m_k_vector() { }

    template< class System >
    inline void do_step( System &system , state_type &x , double t , const double dt )
    {
        fusion::for_each( m_stages , calculate_stage< System >( system , x , m_x_tmp , m_k_vector , t , dt ) );
    }

private:

    const coef_a_type m_a;
    const coef_b_type m_b;
    const coef_c_type m_c;
    const stage_vector m_stages;
    state_type m_x_tmp;
    state_type m_k_vector[num_of_stages];
};

Die einzelnen Stages werden hier auch wieder in einer eigenen Struktur gespeichert: stage_vector . Der Funktor calculate_stage übernimmt die Berechnung einer einzelnen Stage.

Performance

Die Performance der beiden Algorithmen ist etwa äquivalent. Außerdem stehen beide Algorithmen den ausgeschriebenen Varianten in nichts hinterher, in einigen Fällen war der Fusion-Algorithmus sogar schneller als der ausgeschriebene. Genaue Werte sind für verschiedene Compiler unterschiedlich, es gibt sogar Unterschiede zwischen gcc 4.4 und gcc 4.5.

Diskussion

Die Vorteile der beiden Ansätze liegen sofort auf der Hand. Man hat nur genau einen Algorithmus zu warten. Damit können praktisch alle klassischen expliziten Runge-Kutta-Solver (wie Euler, RK4, Midpoint, Dormand-Prince 5-4, Dormand Prince 8-5-3, Fehlberg, ...) geschrieben werden. Außerdem stehen diese Algorithmen den ausgeschriebenen Algorithmen in Performance in nichts nach. Man kann dieses Konzept natürlich auch einfach auf die Stepper mit Fehlerberechnung oder auf Stepper mit kontinuierlichem Output erweitern.

Die Nachteile bei dieser Art von Programmierung sind gering. Der Code ist schwerer zu lesen als klassischer C/C++-Code und es wird auch relativ viel Boilerplate erzeugt. Allerdings wird die dadurch längere Entwicklungszeit durch die Allgemeinheit der Verfahren mehr als wettgemacht.

Selbstverständlich erheben die beiden Versionen des Templatemetaprogrammierungs-Runge-Kutta-Verfahrens keinen Anspruch auf Allgemeingültigkeit. Sie wurden im Rahmen von Boost.Odeint entwickelt und werden dort zum Einsatz kommen. Die Erweiterungsmöglichkeiten dazu sind vielzählig:

Abstrahierung des Zustandstyps. In den obigen beiden Beispielen wird der Zustand der Differentialgleichung durch std::tr1::array< double , N > beschrieben. Das kann selbstverständlich verallgemeinert werden.
Benutzungen von beliebigen Floating-Point-Typen wie float , std::complex< double > , GMP, ...
Erweiterung auf Solver mit Fehlerabschätzung und Dense-Output für Schrittweitensteuerung.
Verallgemeinerung der Algebra.
Höhere Ordnungen der Algebra. Im Moment werden nur Stepper bis Ordnung 4 unterstützt. Das ist trivial erweiterbar auf höhere Ordnungen. Auch können diese Erweiterungen durch Präprozessormakros automatisch erzeugt werden.

Über Kritik, Anregungen, Hinweise sowie Verbesserungsvorschläge würde ich mich sehr freuen.

Literatur

David Abrahams, Aleksey Gurtovoy, C++ Template Metaprogramming: Concepts, Tools, and Techniques from Boost and Beyond (Addison-Wesley, 2004).

Andrei Alexandrescu, Modern C++ Design, Generic Programming and Design Patterns Applied (Addison-Wesley, 2001)

Press William H et al., Numerical Recipes 3rd Edition: The Art of Scientific Computing, 3rd ed. (Cambridge University Press, 2007).

Ernst Hairer, Syvert P. Nørsett, and Gerhard Wanner, Solving Ordinary Differential Equations I: Nonstiff Problems, 2nd ed. (Springer, Berlin, 2009).

Milton Abramowitz, Irene A. Stegun, Handbook of Mathematical Functions (New York: Dover Publications, 1972). (online verfügbar)

headmyshoulder

Super, vielen Dank. Ich hab leider gerade sehr viel um die Ohren. Ich hoffe, dass ich dass morgen das Mergen schaffe.

headmyshoulder

So, ich hab alle Korrekturen eingebaut. Nochmals vielen Dank an die Korrekturleser and die Kommentare. Damit kann der Artikel jetzt veröffentlich werden, oder fehlt noch etwas?

GPC

Ne, passt alles. Du musst die Änderungen jetzt eben noch in den Artikel mit dem [E] Tag einpflegen Weil das ist ja der, den ich veröffentliche.
Dann noch eine veröffentlichbare Version von deinem Autorenprofil (extra Thread) und wir sind im Geschäft

C++ Forumbot

Dieser Thread wurde von Moderator/in GPC aus dem Forum Die Redaktion in das Forum Archiv verschoben.

Im Zweifelsfall bitte auch folgende Hinweise beachten:
C/C++ Forum :: FAQ - Sonstiges :: Wohin mit meiner Frage?

Dieses Posting wurde automatisch erzeugt.