La méthode d'Euler explicite est trop lente pour un problème de réaction-diffusion

10

Je résous le système de réaction-diffusion de Turing avec le code C ++ suivant. C'est trop lent: pour une texture de 128x128 pixels, le nombre d'itérations acceptable est de 200 - ce qui entraîne un délai de 2,5 secondes. J'ai besoin de 400 itérations pour obtenir une image intéressante - mais 5 secondes d'attente, c'est trop. De plus, la taille de la texture devrait être en fait 512x512 - mais cela se traduit par un temps d'attente énorme. Les appareils sont iPad, iPod.

Y a-t-il une chance de le faire plus rapidement? La méthode Euler converge lentement (wikipedia) - une méthode plus rapide permettrait de réduire le nombre d'itérations?

EDIT: Comme l'a souligné Thomas Klimpel, les lignes: "if (m_An [i] [j] <0.0) {...}", "if (m_Bn [i] [j] <0.0) {...}" retardent la convergence: après suppression, une image significative apparaît après 75 itérations . J'ai commenté les lignes dans le code ci-dessous.

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
la source
De plus, je tiens à mentionner qu'il est préférable de ne pas transposer les questions, car il semble que vous ayez posé des questions très similaires ici et ici .
Godric Seer
Avez-vous déjà vu le travail de Greg Turk à ce sujet, par hasard?
JM
@JM: Pas encore. Je viens d'essayer d'exécuter son code: il nécessite un serveur X avec PseudoColor, c'est-à-dire une profondeur de couleur de 8 bits. Je pense que je ne peux pas fournir cela sur OSX. J'ai essayé différents serveurs VNC mais pas de chance.
AllCoder du
Je pense que vous devriez toujours être en mesure d'adapter l'approche de Turk à l'affaire en question; les modèles de réaction-diffusion semblent être assez utilisés de nos jours en infographie.
JM
1
Je me trompe peut-être, mais la partie avec m_An [i] [j] = 0,0; pourrait en fait ajouter un élément à ce système qui ne peut pas être modélisé par une équation différentielle avec un côté droit continu. Cela rend un peu difficile de trouver un solveur plus rapide.
Thomas Klimpel

Réponses:

9

Vous semblez limité par la stabilité, ce qui est attendu car la diffusion est raide lorsque vous affinez la grille. Les bonnes méthodes pour les systèmes rigides sont au moins partiellement implicites. Cela prendra un certain effort, mais vous pouvez implémenter un algorithme multigrille simple (ou utiliser une bibliothèque) pour résoudre ce système avec un coût inférieur à dix "unités de travail" (essentiellement le coût d'un de vos pas de temps). Lorsque vous affinez la grille, le nombre d'itérations n'augmentera pas.

Jed Brown
la source
Si seule la diffusion était rigide ici, il pourrait utiliser une méthode ADI comme Douglas-Gunn et tout irait bien. Cependant, dans ma propre expérience, la partie réactionnelle est souvent bien pire en ce qui concerne la rigidité en plus d'être mal non linéaire.
Thomas Klimpel
1
ADI a malheureusement une localité de mémoire terrible. Notez également que la réaction peut être traitée implicitement, que la diffusion soit ou non. Sous le raffinement de la grille, la diffusion finira par devenir dominante, mais nous ne pouvons pas dire où est le seuil sans connaître les constantes.
Jed Brown
Un exemple de code implémentant Euler à l'envers pour cela (en Python) est ici: scicomp.stackexchange.com/a/2247/123
David Ketcheson
@DavidKetcheson: L'utilisation de méthodes implicites nécessite de résoudre une équation? C'est pourquoi il y a linalg.spsolve () dans le code?
AllCoder du
1
@AllCoder Oui, cela nécessite une résolution, mais la résolution peut être effectuée beaucoup plus rapidement que toutes les étapes de temps requises pour qu'une méthode explicite soit stable.
Jed Brown
2

D'un point de vue pratique: le processeur A5 n'est pas très puissant, vous pouvez donc attendre quelques itérations matérielles, ou si votre ipod / ipad va être connecté à internet, résoudre votre problème à distance ou dans le cloud.

fcruz
la source
Je suis surpris de voir à quel point l'A5 offre peu de puissance. Comment Pages, Safari et autres grandes applications peuvent-elles si bien fonctionner? J'ai besoin de générer des images abstraites aléatoires, j'ai pensé que la morphogenèse serait assez simple ..
AllCoder
Eh bien, A5 est un processeur économe en énergie optimisé pour le Web et la vidéo (Pages, Safari, etc.). En revanche, la plupart des charges de travail numériques effectuent des tonnes d'opérations en virgule flottante et de mouvements de données, ces fonctionnalités ne sont pas au centre des processeurs mobiles à faible puissance.
fcruz
0

Euler converge lentement par rapport à d'autres méthodes, mais je ne pense pas que ce soit ce qui vous intéresse. Si vous recherchez simplement des images "intéressantes", augmentez la taille de votre pas de temps et prenez moins d'itérations. Le problème, comme le souligne Jed, est que la méthode euler explicite a des problèmes de stabilité avec de grands pas de temps par rapport à la taille de la grille. plus votre grille est petite (c'est-à-dire plus la résolution de votre image est élevée), plus votre pas de temps doit être petit pour en tenir compte.

Par exemple, en utilisant euler implicite au lieu d'explicite, vous ne gagnez aucun ordre de convergence, mais la solution aura une stabilité inconditionnelle, permettant des pas de temps beaucoup plus grands. Les méthodes implicites sont plus compliquées à implémenter et prennent plus de calcul par pas de temps, mais vous devriez voir des gains bien au-delà en prenant moins d'étapes au total.

Godric Seer
la source
Ce problème est limité par la stabilité, donc augmenter simplement la taille du pas de temps ne fonctionnera pas.
Jed Brown
Si je change 0,01 par exemple 0,015, alors j'obtiens une "concentration de chimie sp. Proche de zéro" à tous les points - c'est-à-dire un carré gris. Voici l'origine de mon code: drdobbs.com/article/print?articleId=184410024
AllCoder
Oui, ce serait le résultat des problèmes de stabilité mentionnés par Jed. Comme il le mentionne dans sa réponse, l'utilisation d'une méthode implicite caractérisée par de meilleures performances de stabilité résoudra ce problème pour vous. Je mettrai à jour ma réponse pour supprimer les informations non pertinentes.
Godric Seer