Metody rozwiązywania nieliniowych układów doradczo-dyfuzyjnych poza Newton-Raphson?

Pracuję nad projektem, w którym mam dwie domeny sprzężone za pomocą adv-diff za pomocą odpowiednich terminów źródłowych (jedna domena dodaje masę, druga odejmuje masę). Dla zwięzłości modeluję je w stanie ustalonym. Równania to standardowe równanie transportowo-dyfuzyjne z terminem źródłowym wyglądającym tak:

$$\frac{\partial c_1}{\partial t} = 0 = \mathcal{F}_1 + \mathcal{Q}_1(c_1,c_2) \\ \frac{\partial c_2}{\partial t} = 0 = \mathcal{F}_2 + \mathcal{Q}_2(c_1,c_2)$$

Gdzie $\mathcal{F}_i$ jest dyfuzyjnym i korzystnym strumieniem dla gatunków $i$, i $\mathcal{Q}_i$ jest terminem źródłowym dotyczącym gatunku $i$.

Udało mi się napisać solver dla mojego problemu za pomocą metody Newtona-Raphsona i całkowicie sparowałem dwie domeny za pomocą macierzy mas blokowych, tj .:

$$F_{coupled} = \left[\begin{array}{c c} A_1 & 0 \\ 0 & A_2 \\ \end{array}\right]\underbrace{ \left[\begin{array}{c} c_{1,i} \\ c_{2,i} \\ \end{array}\right] }_{x_i} - \left[\begin{array}{c} b_1(c_{1,i}, c_{2,i}) \\ b_2(c_{1,i}, c_{2,i}) \\ \end{array}\right]$$

Termin $F_{coupled}$ służy do określania macierzy jakobiańskiej i aktualizacji obu $c_1$ i $c_2$:

$$\mathcal{J}(x_i) \left[x_{i+1} - x_i \right] = -\mathcal{F}_{coupled}$$

lub

$$x_{i+1} = x_i - \left(\mathcal{J}(x_i)\right)^{-1} \mathcal{F}_{coupled}$$

Aby przyspieszyć, nie obliczam jakobianu przy każdej iteracji - w tej chwili bawię się co pięć iteracji, które wydają się działać wystarczająco dobrze i utrzymywać rozwiązanie na stałym poziomie.

Problem polega na tym: przeniosę się do większego systemu, w którym obie domeny są w 2D / 2.5D, a obliczenie macierzy jakobińskiej szybko wyczerpie moje dostępne zasoby komputerowe. Buduję ten model do późniejszego wykorzystania w ustawieniach optymalizacyjnych, więc nie mogę też być za kierownicą przy każdej iteracji, dostosowując współczynnik tłumienia itp.

Czy mam rację, szukając gdzie indziej bardziej niezawodnego i algorytmu dla mojego problemu, czy jest to tak dobre, jak to możliwe? Zajrzałem trochę do quasi-linearyzacji, ale nie jestem pewien, jak ma zastosowanie w moim systemie.

Czy są jakieś inne zręczne algorytmy, które mogłem przeoczyć, które mogą rozwiązać układ równań nieliniowych bez uciekania się do ponownego obliczania jakobianów jako przestępstwa?

Zakładam, że ograniczeniem w 2D i 3D jest przechowywanie jakobianów.

Jedną z opcji jest zachowanie pochodnych czasu i użycie wyraźnego „pseudo” krokowego czasu, aby przejść do stanu ustalonego. Zwykle liczba CFL potrzebna do systemów dyfuzyjnych i reaktywnych może być zbyt mała. Możesz wypróbować nieliniową multigrid (jest to również nazywane multigrid z pamięcią o pełnej aproksymacji) i lokalnym krokiem w celu przyspieszenia konwergencji.

Inną opcją jest użycie w pełni niejawnego schematu, tak jak teraz, ale nie przechowywanie globalnego jakobianu. Możesz użyć niejawnego schematu wolnego od macierzy.

$$DF(u^n)\, \delta u^n = -F(u^n)$$ (gdzie $DF$ jest jakobianem) można rozwiązać za pomocą solvera do podprzestrzeni Kryłowa, takiego jak GMRES i BiCGStab, wykorzystując fakt, że $$DF(u^n)\, \delta u \approx \frac{F(u^n+\epsilon\frac{\delta u}{\Vert\delta u\Vert}) - F(u^n)}{\epsilon}.$$ Wynika to z faktu, że GMRES i BiCGStab nie wymagają matrycy LHS $A$, muszą tylko móc obliczyć swój produkt $Ax$ dany wektor $x$.

Teraz z odpowiednią wartością $\epsilon$ (zwykle około $10^{-7}$dla pływaków podwójnej precyzji) możesz wykonać pętlę Newtona bez obliczania lub przechowywania jakobianów. Wiem na pewno, że ta technika jest używana do rozwiązywania niektórych nietrywialnych przypadków w obliczeniowej dynamice płynów. Należy jednak pamiętać, że liczba ocen funkcji$F$ będzie więcej niż w technice przechowywania macierzy, zamiast wymagać produktu macierz-wektor.

Inną rzeczą do odnotowania jest to, że jeśli twój system jest taki, że potrzebny jest potężny warunek wstępny (tj. Jacobi lub blok-Jacobi nie wystarczą), możesz spróbować użyć wyżej wspomnianej metody jako płynniejszej w schemacie wielosieciowym. Jeśli chcesz wypróbować punktowy lub blokowy warunek wstępny Jacobiego, możesz obliczyć i zapisać tylko ukośne elementy lub ukośne bloki jakobianu, co nie jest wiele. Chciałbym również wspomnieć, że wdrożenie warunku Gaussa-Seidla lub SSOR może być możliwe bez wyraźnego przechowywania jakobianu. W tym artykule opisano implementację wolnego od macierzy GMRES wstępnie wyposażonego w pozbawiony macierzy symetryczny Gaussa-Seidela w kontekście obliczeniowej dynamiki płynów.

Z mojego doświadczenia z równaniami Naviera-Stokesa można sobie radzić bardzo dobrze bez całkowicie ukrytych schematów.

Jeśli potrzebujesz tylko szybkiego schematu numerycznego dla rozwiązania ewolucji czasu, spójrz na schematy IMEX (niejawne); patrz np. ten artykuł Aschera, Ruuth, Spiteri Implicit-Explicit Metody Runge-Kutta dla zależnych od czasu równań częściowych różniczkowych .

Możesz nawet spróbować po prostu użyć jawnego schematu integracji czasowej wysokiego poziomu z kontrolą wielkości kroku (jak Matlab ODE45). Możesz jednak wpaść w kłopoty z powodu sztywności twojego systemu, która pochodzi z części dyfuzyjnej. Na szczęście część dyfuzyjna jest liniowa, dzięki czemu można ją leczyć niejawnie (taka jest idea schematów IMEX).

Na początku zastanawiałem się nad dodaniem tylko uwagi, ale miejsca było za mało, dlatego dodam krótki opis moich doświadczeń z tym tematem.

Po pierwsze, spójrz na swoją notację $F_{coupled}$ Przypuszczam, że nie widzę sprzężonej formy $b_1$ i $b_2$ będą zależne od $c_{1,i}$ i $c_{2,i}$. Ponadto jeśli$A_1$ i $A_2$ są reprezentacjami macierzowymi przybliżeń $\mathcal{F}_1$ i $\mathcal{F}_2$ to nie powinny zależeć tylko od $c_i$, ale także wartości sąsiada, ale może to być tylko błędne zrozumienie notacji.

Jako komentarz ogólny chciałbym dodać, że użycie analitycznego jakobianu wydaje się być jedynym sposobem na uzyskanie kwadratowej zbieżności nieliniowego solwera iteracyjnego (tzn. W twoim przypadku solvera Newtona-Raphsona). Czy zaobserwowałeś to w swoim przypadku? Jest to dość ważne, ponieważ w przeciwnym razie może dojść do nieporozumień w przybliżeniach (linearyzacja).

We wszystkich aplikacjach, w które byłem zaangażowany (niektóre z nich obejmowały obliczenia na dużą skalę) nigdy nie mieliśmy problemów z czasochłonnością składania jakobianów, najbardziej czasochłonnym problemem zawsze było zastosowanie liniowego solvera. Analityczny jakobian (jeśli jest dostępny) zawsze był w aplikacjach, nad którymi pracowałem nad preferowanym wyborem ze względu na kwadratową zbieżność. W nielicznych przypadkach taki nieliniowy solver wytwarza macierz, która powoduje problemy z konwergencją iteracyjnego solwera liniowego, więc próbowaliśmy zastosować prostszą linearyzację niż analityczna metoda jakobowska, aby pomóc liniowemu solverowi. Taki kompromis między zachowaniem nieliniowego i liniowego rozwiązania algebraicznego w zależności od linearyzacji nieliniowego układu algebraicznego zawsze był trudny i nie mogłem dać ogólnej rekomendacji.

Ale masz rację, że wadą (lub właściwością) analitycznego jakobianu dla systemu PDE jest to, że wytwarza on sprzężony system algebraiczny, więc jeśli oddzielisz taki system w jakikolwiek sposób, np. Rozwiązując osobno przybliżenie każdego PDE, powiedzmy, iteracyjny podział następnie tracisz kwadratową zbieżność globalnego solvera. Ale przynajmniej, jeśli rozwiążesz osobno każdy dyskretny (oddzielony) PDE, możesz ponownie przyspieszyć rozwiązanie tego konkretnego problemu, stosując metodę Newtona-Raphsona.