Explizite Euler-Methode zu langsam für Reaktionsdiffusionsprobleme

10

Ich löse Turings Reaktionsdiffusionssystem mit folgendem C ++ - Code. Es ist zu langsam: Für eine Textur mit 128 x 128 Pixeln beträgt die akzeptable Anzahl von Iterationen 200 - was zu einer Verzögerung von 2,5 Sekunden führt. Ich benötige 400 Iterationen, um ein interessantes Bild zu erhalten - aber 5 Sekunden Warten sind zu viel. Außerdem sollte die Größe der Textur tatsächlich 512 x 512 betragen - dies führt jedoch zu einer enormen Wartezeit. Die Geräte sind iPad, iPod.

Gibt es eine Chance, dies schneller zu tun? Die Euler-Methode konvergiert langsam (Wikipedia) - eine schnellere Methode würde es ermöglichen, die Anzahl der Iterationen zu verringern?

EDIT: Wie Thomas Klimpel betonte, lauten die Zeilen: "if (m_An [i] [j] <0.0) {...}", "if (m_Bn [i] [j] <0.0) {...}" verzögern die Konvergenz: Nach dem Entfernen erscheint nach 75 Iterationen ein aussagekräftiges Bild . Ich habe die Zeilen im Code unten auskommentiert.

void TuringSystem::solve( int iterations, double CA, double CB ) {
    m_iterations = iterations;
    m_CA = CA;
    m_CB = CB;

    solveProcess();
}

void set_torus( int & x_plus1, int & x_minus1, int x, int size ) {
    // Wrap "edges"
    x_plus1 = x+1;
    x_minus1 = x-1;
    if( x == size - 1 ) { x_plus1 = 0; }
    if( x == 0 ) { x_minus1 = size - 1; }
}

void TuringSystem::solveProcess() {
    int n, i, j, i_add1, i_sub1, j_add1, j_sub1;
    double DiA, ReA, DiB, ReB;

    // uses Euler's method to solve the diff eqns
    for( n=0; n < m_iterations; ++n ) {
        for( i=0; i < m_height; ++i ) {
            set_torus(i_add1, i_sub1, i, m_height);

            for( j=0; j < m_width; ++j ) {
                set_torus(j_add1, j_sub1, j, m_width);

                // Component A
                DiA = m_CA * ( m_Ao[i_add1][j] - 2.0 * m_Ao[i][j] + m_Ao[i_sub1][j]   +   m_Ao[i][j_add1] - 2.0 * m_Ao[i][j] + m_Ao[i][j_sub1] );
                ReA = m_Ao[i][j] * m_Bo[i][j] - m_Ao[i][j] - 12.0;
                m_An[i][j] = m_Ao[i][j] + 0.01 * (ReA + DiA);
                // if( m_An[i][j] < 0.0 ) { m_An[i][j] = 0.0; }

                // Component B
                DiB = m_CB * ( m_Bo[i_add1][j] - 2.0 * m_Bo[i][j] + m_Bo[i_sub1][j]   +   m_Bo[i][j_add1] - 2.0 * m_Bo[i][j] + m_Bo[i][j_sub1] );
                ReB = 16.0 - m_Ao[i][j] * m_Bo[i][j];
                m_Bn[i][j] = m_Bo[i][j] + 0.01 * (ReB + DiB);
                // if( m_Bn[i][j] < 0.0 ) { m_Bn[i][j]=0.0; }
            }
        }

        // Swap Ao for An, Bo for Bn
        swapBuffers();
    }
}
AllCoder
quelle
Außerdem möchte ich erwähnen, dass es bevorzugt wird, dass Sie keine Fragen über Kreuz stellen, da Sie anscheinend sowohl hier als auch hier sehr ähnliche Fragen gestellt haben .
Godric Seer
Haben Sie vielleicht schon Greg Turk's Arbeit dazu gesehen?
JM
@JM: Noch nicht. Ich habe gerade versucht, seinen Code auszuführen: Er erfordert einen X-Server mit PseudoColor, dh eine 8-Bit-Farbtiefe. Ich glaube, ich kann dies unter OSX nicht bereitstellen. Ich habe verschiedene VNC-Server ausprobiert, aber kein Glück.
AllCoder
Ich denke, Sie sollten immer noch in der Lage sein, den Ansatz von Turk an die vorliegende Angelegenheit anzupassen. Reaktionsdiffusionsmuster scheinen heutzutage in der Computergrafik ein gutes Stück verwendet zu werden.
JM
1
Ich könnte mich irren, aber der Teil mit m_An [i] [j] = 0.0; könnte diesem System tatsächlich ein Element hinzufügen, das nicht durch eine Differentialgleichung mit einer kontinuierlichen rechten Seite modelliert werden kann. Dies macht es etwas schwierig, einen schnelleren Löser zu finden.
Thomas Klimpel

Antworten:

9

Sie scheinen durch die Stabilität eingeschränkt zu sein, was zu erwarten ist, da die Diffusion beim Verfeinern des Gitters steif ist. Gute Methoden für steife Systeme sind zumindest teilweise implizit. Es wird einige Mühe kosten, aber Sie können einen einfachen Multigrid-Algorithmus implementieren (oder eine Bibliothek verwenden), um dieses System mit Kosten von weniger als zehn "Arbeitseinheiten" (im Wesentlichen den Kosten eines Ihrer Zeitschritte) zu lösen. Wenn Sie das Raster verfeinern, wird die Anzahl der Iterationen nicht erhöht.

Jed Brown
quelle
Wenn hier nur die Diffusion steif wäre, könnte er eine ADI-Methode wie Douglas-Gunn verwenden und alles wäre in Ordnung. Nach meiner eigenen Erfahrung ist der Reaktionsteil jedoch in Bezug auf die Steifheit oft viel schlechter als schlecht nichtlinear.
Thomas Klimpel
1
ADI hat leider eine schreckliche Speicherlokalität. Beachten Sie auch, dass die Reaktion implizit behandelt werden kann, unabhängig davon, ob es sich um eine Diffusion handelt. Bei der Verfeinerung des Gitters wird die Diffusion schließlich dominant, aber wir können nicht sagen, wo der Schwellenwert liegt, ohne die Konstanten zu kennen.
Jed Brown
Beispielcode für die Implementierung von Backward Euler (in Python) finden Sie hier: scicomp.stackexchange.com/a/2247/123
David Ketcheson
@ DavidKetcheson: Um implizite Methoden zu verwenden, muss eine Gleichung gelöst werden. Aus diesem Grund enthält der Code linalg.spsolve ()?
AllCoder
1
@AllCoder Ja, es ist eine Lösung erforderlich, aber die Lösung kann viel schneller durchgeführt werden als alle Zeitschritte, die erforderlich sind, damit eine explizite Methode stabil ist.
Jed Brown
2

Aus praktischer Sicht: Der A5-Prozessor ist nicht sehr leistungsfähig, sodass Sie einige HW-Iterationen abwarten können. Wenn Ihr iPod / iPad mit dem Internet verbunden werden soll, lösen Sie Ihr Problem remote oder in der Cloud.

fcruz
quelle
Ich bin überrascht, wie wenig Leistung der A5 bietet. Wie können Pages, Safari und andere große Anwendungen so gut funktionieren? Ich muss zufällige, abstrakte Bilder erzeugen, dachte, dass die Morphogenese einfach genug sein wird.
AllCoder
Nun, A5 ist ein energieeffizienter Prozessor, der für Web und Video (Pages, Safari usw.) optimiert ist. Im Gegensatz dazu führen die meisten numerischen Workloads Tonnen von Gleitkommaoperationen und Datenbewegungen aus. Diese Funktionen stehen nicht im Mittelpunkt mobiler Prozessoren mit geringem Stromverbrauch.
Fcruz
0

Euler konvergiert zwar langsam in Bezug auf andere Methoden, aber ich denke nicht, dass Sie daran interessiert sind. Wenn Sie nur nach "interessanten" Bildern suchen, vergrößern Sie Ihren Zeitschritt und führen Sie weniger Iterationen durch. Jed weist darauf hin, dass das Problem darin besteht, dass die explizite Eulermethode Stabilitätsprobleme mit großen Zeitschritten in Bezug auf die Gittergröße aufweist. Je kleiner Ihr Raster ist (dh je höher die Auflösung Ihres Bildes ist), desto kleiner muss Ihr Zeitschritt sein, um dies zu berücksichtigen.

Wenn Sie beispielsweise implizites Euler anstelle von explizit verwenden, erhalten Sie keine Konvergenzreihenfolgen, aber die Lösung weist eine bedingungslose Stabilität auf, die viel größere Zeitschritte ermöglicht. Implizite Methoden sind komplizierter zu implementieren und erfordern mehr Berechnungen pro Zeitschritt. Sie sollten jedoch weit darüber hinaus Gewinne erzielen, wenn Sie insgesamt weniger Schritte ausführen.

Godric Seer
quelle
Dieses Problem ist durch die Stabilität begrenzt, sodass eine einfache Erhöhung der Zeitschrittgröße nicht funktioniert.
Jed Brown
Wenn ich 0,01 auf zB 0,015 ändere, erhalte ich an allen Punkten eine "Konzentration von chem. Sp. Nahe Null" - dh ein graues Quadrat. Hier ist der Ursprung meines Codes: drdobbs.com/article/print?articleId=184410024
AllCoder
Ja, das wäre ein Ergebnis der von Jed erwähnten Stabilitätsprobleme. Wie er in seiner Antwort erwähnt, wird dieses Problem durch die Verwendung einer impliziten Methode, die durch eine bessere Stabilitätsleistung gekennzeichnet ist, für Sie behoben. Ich werde meine Antwort aktualisieren, um die irrelevanten Informationen zu entfernen.
Godric Seer