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

10

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();
    }
}
AllCoder
źródło
Chciałbym również wspomnieć, że lepiej nie zadawać pytań, ponieważ wygląda na to, że zadałeś bardzo podobne pytania zarówno tutaj, jak i tutaj .
Godric Seer,
Czy przypadkiem widziałeś już prace Grega Turka ?
JM
@JM: Jeszcze nie. Właśnie próbowałem uruchomić jego kod: wymaga X serwera z PseudoColor, tj. 8-bitowa głębia kolorów. Myślę, że nie mogę tego zapewnić w OSX. Próbowałem różnych serwerów VNC, ale bez powodzenia.
AllCoder,
Myślę, że nadal powinieneś być w stanie dostosować podejście Turk'a do omawianej sprawy; wzory dyfuzji reakcji wydają się być obecnie dość często używane w grafice komputerowej.
JM
1
Mogę się mylić, ale część z m_An [i] [j] = 0,0; może faktycznie dodać element do tego układu, którego nie można modelować równaniem różniczkowym z ciągłą prawą stroną. Utrudnia to wymyślenie szybszego solvera.
Thomas Klimpel

Odpowiedzi:

9

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.

Jed Brown
źródło
Jeśli tylko dyfuzja byłaby tutaj sztywna, mógłby użyć metody ADI, takiej jak Douglas-Gunn i wszystko byłoby w porządku. Jednak z własnego doświadczenia wynika, że ​​część reakcji jest często znacznie gorsza pod względem sztywności, a ponadto jest źle nieliniowa.
Thomas Klimpel
1
Niestety ADI ma okropną lokalizację pamięci. Należy również pamiętać, że reakcję można leczyć niejawnie, niezależnie od tego, czy jest to dyfuzja. W ramach udoskonalenia siatki dyfuzja w końcu stanie się dominująca, ale nie możemy ustalić, gdzie jest próg, nie znając stałych.
Jed Brown
Przykładowy kod implementujący do tyłu Euler (w Pythonie) znajduje się tutaj: scicomp.stackexchange.com/a/2247/123
David Ketcheson
@DavidKetcheson: Korzystanie z metod niejawnych wymaga rozwiązania równania? Dlatego w kodzie jest linalg.spsolve ()?
AllCoder
1
@AllCoder Tak, wymaga rozwiązania, ale rozwiązanie można wykonać znacznie szybciej niż wszystkie kroki czasowe wymagane dla stabilnej metody jawnej.
Jed Brown
2

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.

fcruz
źródło
Dziwi mnie, jak mała moc oferuje A5. W jaki sposób Pages, Safari i inne duże aplikacje działają tak dobrze? Muszę generować losowe, abstrakcyjne obrazy, myślałem, że morfogeneza będzie wystarczająco prosta ..
AllCoder
A5 to energooszczędny procesor zoptymalizowany pod kątem Internetu i wideo (Pages, Safari itp.). W przeciwieństwie do tego, większość obciążeń numerycznych wykonuje mnóstwo operacji zmiennoprzecinkowych i przenoszenia danych, te funkcje nie są przedmiotem zainteresowania mobilnych procesorów o niskiej mocy.
fcruz
0

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.

Godric Seer
źródło
Ten problem jest ograniczony przez stabilność, więc po prostu zwiększenie rozmiaru kroku czasu nie będzie działać.
Jed Brown
Jeśli zmienię 0,01 na np. 0,015, to otrzymam „stężenie chem. Sp. Blisko zera” we wszystkich punktach - tzn. Szary kwadrat. Oto pochodzenie mojego kodu: drdobbs.com/article/print?articleId=184410024
AllCoder
Tak, byłby to wynik problemów ze stabilnością, o których wspominał Jed. Jak wspomina w swojej odpowiedzi, zastosowanie metody niejawnej charakteryzującej się lepszą wydajnością stabilności rozwiąże ten problem. Zaktualizuję swoją odpowiedź, aby usunąć nieistotne informacje.
Godric Seer,