Jawna metoda Eulera jest zbyt wolna, aby rozwiązać problem dyfuzji reakcji

Rozwiązuję układ dyfuzyjny reakcji Turinga za pomocą następującego kodu C ++. Jest zbyt wolny: dla tekstury 128 x 128 pikseli dopuszczalna liczba iteracji wynosi 200 - co powoduje 2,5 sekundy opóźnienia. Potrzebuję 400 iteracji, aby uzyskać interesujący obraz - ale 5 sekund oczekiwania to za dużo. Ponadto rozmiar tekstury powinien w rzeczywistości wynosić 512 x 512 - ale powoduje to długi czas oczekiwania. Urządzeniami są iPad, iPod.

Czy jest jakaś szansa, aby zrobić to szybciej? Metoda Eulera zbiega się powoli (wikipedia) - czy szybsza metoda pozwoliłaby na zmniejszenie liczby iteracji?

EDYCJA: Jak zauważył Thomas Klimpel, linie: „if (m_An [i] [j] <0,0) {...}”, „if (m_Bn [i] [j] <0,0) {...}” opóźniają konwergencję: po usunięciu znaczący obraz pojawia się po 75 iteracjach . Skomentowałem poniższe wiersze w kodzie.

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();
    }
}

Wygląda na to, że jesteś ograniczony przez stabilność, której oczekuje się, ponieważ dyfuzja jest sztywna podczas udoskonalania siatki. Dobre metody dla sztywnych układów są przynajmniej częściowo niejawne. To zajmie trochę wysiłku, ale możesz wdrożyć prosty algorytm wielosiatkowy (lub użyć biblioteki), aby rozwiązać ten system kosztem mniejszym niż dziesięć „jednostek roboczych” (zasadniczo koszt jednego z twoich kroków czasowych). Po poprawieniu siatki liczba iteracji nie wzrośnie.

Z praktycznego punktu widzenia: procesor A5 nie jest tak potężny, więc możesz poczekać kilka iteracji HW lub jeśli twój ipod / ipad będzie podłączony do Internetu, rozwiąż swój problem zdalnie lub w chmurze.

Euler zbiega się powoli w stosunku do innych metod, jednak nie sądzę, że jesteś tym zainteresowany. Jeśli szukasz tylko „interesujących” zdjęć, zwiększ rozmiar swojego kroku czasowego i wykonuj mniej iteracji. Problem, jak zauważa Jed, polega na tym, że jawna metoda eulera ma problemy ze stabilnością z dużymi krokami czasowymi w stosunku do rozmiaru siatki. im mniejsza jest twoja siatka (tj. wyższa rozdzielczość twojego obrazu), tym mniejszy musi być twój krok czasowy.

Na przykład, używając domyślnego eulera zamiast jawnego, nie zyskujesz żadnych rzędów zbieżności, ale rozwiązanie będzie miało bezwarunkową stabilność, umożliwiając znacznie większe kroki czasowe. Metody niejawne są bardziej skomplikowane w implementacji i wymagają więcej obliczeń na krok, ale powinieneś zobaczyć znacznie więcej, biorąc pod uwagę mniej kroków ogółem.