Czy można rozwiązać nieliniowe PDE bez iteracji Newtona-Raphsona?

Próbuję zrozumieć niektóre wyniki i doceniłbym kilka ogólnych uwag na temat rozwiązywania problemów nieliniowych.

Równanie Fishera (nieliniowe PDE z dyfuzją reakcyjną),

u_{t} = re u_{x x} + β u (1 - u) = fa (u)

$u_t = du_{xx} + \beta u (1 - u) = F(u)$

w formie dyskretnej,

u_{jot}^{'} = L. u + β u_{jot} (1 - u_{jot}) = fa (u)

$u_j^{\prime} = \boldsymbol{L}\boldsymbol{u} + \beta u_j (1 - u_j) = F(\boldsymbol{u})$

gdzie $\boldsymbol{L}$ jest operatorem różnicowym, a $\boldsymbol{u}=(u_{j-1}, u_j, u_{j+1})$ to szablon dyskretyzacji.

metoda

Chcę zastosować system domyślny, ponieważ potrzebuję stabilności i nieograniczonego czasu. W tym celu używam metody $\theta$ , (zauważ, że $\theta=1$ daje w pełni niejawny schemat, a $\theta=0.5$ daje schemat trapezoidalny lub „Crank-Nicolson”),

u_{j}^{'} = θ F (u^{n + 1}) + (1 - θ) F (u^{n})

$u_{j}^{\prime} = \theta F(\boldsymbol{u}^{n+1}) + (1-\theta) F(\boldsymbol{u}^{n})$

Jednak w przypadku problemów nieliniowych nie można tego zrobić, ponieważ równania nie można zapisać w postaci liniowej.

Aby obejść ten problem, badałem dwa podejścia numeryczne,

Metoda IMEX

$u_{jot}^{'} = \underset{θ - metoda dyfuzji}{\underset{⏟}{θ L. u^{n + 1} + (1 - θ) L. u^{n}}} + \underset{W pełni wyraźny termin reakcji}{\underset{⏟}{β u_{jot}^{n} (1 - u_{jot}^{n})}}$ $u_j^{\prime} = \underbrace{\theta\boldsymbol{L}\boldsymbol{u^{n+1}} + (1-\theta)\boldsymbol{L}\boldsymbol{u^{n}}}_{\theta-\text{method diffusion term}} + \underbrace{\beta u_j^{n} (1 - u_j^{n})}_{\text{Fully explicit reaction term}}$
Najbardziej oczywistą drogą jest zignorowanie nieliniowej części terminu reakcji i po prostu zaktualizowanie terminu reakcji z najlepszą możliwą wartością, tj. Tą z poprzedniego etapu czasowego. Wynikiem tego jest metoda IMEX.
Solver Newtona

ν^{k + 1} = ν^{k} - (I - θ τ A^{n})^{- 1} (ν^{k} - u_{n} - (1 - θ) τ F (w^{n}) - θ τ F (w^{n + 1}))

$\nu^{k+1} = \nu^{k} - (I - \theta\tau A^{n})^{-1} \left( \nu^{k} - u_{n} - (1-\theta) \tau F(w^n) - \theta\tau F(w^{n+1}) \right)$

Równanie pełnej -method można rozwiązać za pomocą iteracji Newtona-Raphsona, aby znaleźć zmienną przyszłego rozwiązania. Gdzie jest indeksem iteracji ( ), a jest jakobską macierzą . W tym przypadku używam symboli dla zmiennych iteracyjnych, tak że różnią się one od rozwiązania równania w punkcie czasu rzeczywistego . W rzeczywistości jest to zmodyfikowany solver Newtona, ponieważ jakobian nie jest aktualizowany przy każdej iteracji. $\theta$ $k$ $k\geq0$ $A^{n}$ $F(w^n)$ $\nu^{k}$ $u^n$

Wyniki

Porównanie równań Fishera metod numerycznych.

Powyższe wyniki są obliczane dla stosunkowo dużego kroku czasowego i pokazują różnicę między podejściem krokowym a pełnym rozwiązaniem iteracyjnym Newtona.

Czego nie rozumiem:

Dziwi mnie, że metoda skokowa czasu działa „OK”, ale z czasem pozostaje w tyle za rozwiązaniem analitycznym. ( Uwaga: jeśli wybrałem mniejszy krok czasowy, wówczas podejście krokowe daje wyniki zamknięte dla modelu analitycznego). Dlaczego podejście krokowe daje rozsądne wyniki równania nieliniowego?
Model Newtona ma się znacznie lepiej, ale z czasem zaczyna przewodzić modelowi analitycznemu. Dlaczego dokładność podejścia Newtona zmniejsza się z czasem? Czy można poprawić dokładność?
Dlaczego istnieje ogólna cecha, że po wielu iteracjach model numeryczny i model analityczny zaczynają się różnić? Czy to tylko dlatego, że przedział czasowy jest zbyt duży, czy zawsze tak się stanie?

— boyfarrell
źródło

Polecam przeczytanie podstawowej analizy błędów solverów ODE, na przykład w Hairer / Nørsett / Wanner, a także niektórych analiz stabilności. Odpowiedzi na większość twoich pytań.

— Guido Kanschat

@boyfarrell, aby uniknąć nieporozumień wśród innych czytelników, powinieneś położyć terminologię w miejscu, w którym wyjaśniasz swoją metodę: 1. IMEX - wyraźny w nieliniowości i ukryty w części liniowej. 2. jest to średnia -schemat, które zazwyczaj wymagają metodę Newtona do rozwiązania dla aktualizacji

θ

$\theta$

— Jan

Cześć @ Jan Myślę, że mam wszystko. Jeszcze raz dziękuję za pomoc.

— boyfarrell

Odpowiedzi:

Zakładam, że przeprowadziłeś dyskretyzację przestrzeni, więc chodzi o rozwiązanie ODE za pomocą schematu numerycznego który przesuwa przybliżenie w chwili obecnej do następnej wartości przy .

{\dot{u}}_{h} (t) = {fa}_{h} (t, u_{h} (t)), na [0, T], u_{h} (0) = α .

$\dot u_h(t) = F_h(t,u_h(t)), \text{ on [0,T] }, u_h(0) = \alpha.$

Φ

$\Phi$

u_{h}^{n}

$u_h^n$

t = t^{n}

$t=t^n$

u_{h}^{n + 1}

$u_h^{n+1}$

t = t^{n + 1} := t^{n} + τ

$t=t^{n+1}:=t^n+\tau$

Następnie twoje pytania odnoszą się do właściwości jawnych , gdzie aktualizacja zapisuje jako

u_{h}^{n + 1} = u_{h}^{n} + Φ_{mi} (t^{n}, τ, u_{h}^{n}),

$u_h^{n+1} = u_h^n + \Phi_e(t^n,\tau,u_h^n),\quad \quad \quad$

niejawne , zapisane jak

u_{h}^{n + 1} = u_{h}^{n} + Φ_{ja} (t^{n}, τ, u_{h}^{n + 1}, u_{h}^{n}), (*)

$u_h^{n+1} = u_h^n + \Phi_i(t^n,\tau,u_h^{n+1},u_h^n), \quad \quad \quad (*)$

lub kombinacja obu („ IMEX ”, patrz odpowiedź @Jeda Browna) jednoetapowych schematów krokowych.

W tym ustawieniu metoda Newtona jest po prostu podejściem do rozwiązania potencjalnie nieliniowego systemu wynikającego z . $u_h^{n+1}$ $(*)$

A moje odpowiedzi opierają się na wynikach analizy numerycznej metod jednoetapowych.

W przypadku korzystania ze schematów konwergencji pod względem kolejności konwergencji nie ma ogólnej korzyści z zastosowania schematów niejawnych (patrz. 2.). Jednak dla sztywnych systemów, np. Twojego systemu zawierającego Laplaciana, istnieją ukryte schematy, które są stabilne bez ograniczeń czasowych. Niemniej jednak teoretycznie w przypadku schematu jawnego uzyskuje się lepsze wyniki przy mniejszych czasowych , o ile samo równanie jest stabilne (np. Odnosząc się do twierdzenia Picarda-Lindelofa, jeśli to Lipshitz w drugim argumencie) i czasu -step nie jest zbyt mały. $F_h$
Możesz znaleźć przykłady, w których wyraźne schematy działają lepiej. (Teoretycznie możesz odwrócić czas w swoim przykładzie, zacząć od wartości końcowej i znaleźć ukrytą i jawną zamianę.) Jeśli sprawisz, że błąd Newtona będzie wystarczająco mały, nadal możesz poprawić dokładność, zmniejszając krok czasowy lub wykorzystując czas -stepujące schematy wyższego rzędu.
Stała w oszacowaniu błędu dla błędu globalnego rośnie wykładniczo wraz z długością przedziału czasowego. Zobacz, np. Tutaj, wyraźny schemat Eulera. Dotyczy to każdej metody jednoetapowej. Ponieważ szacunek jest typu , , mniejszy krok jedynie opóźnia ten efekt. $C$ $err \approx C \tau^p$ $p>0$ $\tau$

Jeszcze kilka uwag i ostateczna odpowiedź:

Schematy IMEX mogą być stosowane do niejawnego traktowania tylko części liniowej, co pozwala uniknąć rozwiązań nieliniowych. Zobacz odpowiedź Jeda Browna.
Crank-Nicolson jest metodą jednoetapową. „Wiele” w metodach wieloetapowych odnosi się do zastosowania szeregu wcześniejszych kroków czasowych do zdefiniowania bieżącej aktualizacji. Np. Jak Jest to bardzo różne od metod jednoetapowych, a także dzielonych lub IMEX, w których aktualizacja jest zdefiniowana nie powtarzając się do poprzednich wartości. $u_{h}^{n + 1} = Φ_{m} (t^{n}, τ, u_{h}^{n + 1}, u_{h}^{n}, u_{h}^{n - 1}) .$ $u_h^{n+1} = \Phi_m(t^n,\tau,u_h^{n+1},u_h^n,u_h^{n-1}).$

Więc moja odpowiedź brzmi: tak , możesz rozwiązać nieliniowe PDE bez metody Newtona. Możesz użyć schematów jawnych, schematów „IMEX” lub tak zwanych metod niejawnie liniowych (np. Metod Rosenbrocka). Możesz także zastosować inne podejścia do rozwiązywania układów z takie jak iteracja w punkcie stałym lub, w szczególnych przypadkach, solwery algebraiczne. $(*)$

— Jan
źródło

Tak, zastosowałem standardowy centralny szablon różnicowy do terminu dyfuzji. Nie mogę użyć jawnego schematu (dla prawdziwego problemu, który chcę rozwiązać), ponieważ stabilny krok czasowy jest nierealnie niewielki. Dlatego badam IMEX lub ukryte opcje. Jeśli chodzi o twój trzeci punkt, aby uniknąć kumulacji błędów, muszę zastosować metody wieloetapowe. Czy schemat Crank-Nicolson, którego użyłem powyżej (z solwerem Newtona), jest sklasyfikowany jako metoda wieloetapowa (ma dwa punkty w czasie)? Byłem zaskoczony, że błąd narastał z czasem przy stosowaniu metody solvera Newtona.

— boyfarrell

Crank-Nicolson jest metodą jednoetapową, ponieważ pisze jako . Nie rozumiem też, dlaczego schematy wieloetapowe powinny unikać akumulacji błędów.

u_{h}^{n + 1} = u_{h}^{n} + Φ (t^{n}, τ^{n}, u_{h}^{n}, u_{h}^{n + 1})

$u_h^{n+1}=u_h^n + \Phi(t^n,\tau^n,u_h^n,u_h^{n+1})$

— Jan

OK, dziękuję za wyjaśnienie na temat metody CN. Tak, interesujące jest to, dlaczego metody wieloetapowe wydają się mieć mniejszą akumulację błędów. Rozumiem teraz, że solver Newtona ma narastający błąd, ponieważ jest to metoda jednoetapowa. Nawiasem mówiąc, wiem, że lubisz Python. Zrobiłem wszystkie powyższe przy użyciu scipy, numpy i matplotlib, gist.github.com/danieljfarrell/6353776

— boyfarrell

Mam usunięto link do papieru przez Trefethen et. glin. o wysokiej kolejności integracji IMEX z mojej odpowiedzi, ponieważ istnieją lepsze odniesienia do poznania schematów IMEX.

— Jan

Krótka odpowiedź

Jeśli potrzebujesz tylko dokładności drugiego rzędu i braku wbudowanej oceny błędów, istnieje szansa, że będziesz zadowolony z podziału Strang: pół kroku reakcji, pełny krok dyfuzji, pół kroku reakcji.

Długa odpowiedź

Dyfuzja reakcji, nawet przy reakcji liniowej, słynie z wykazywania błędu podziału. Rzeczywiście, może być znacznie gorzej, w tym „konwergowanie” do niepoprawnych stanów ustalonych, mylenie stanów ustalonych z cyklami granicznymi, mylenie konfiguracji stabilnych i niestabilnych i wiele innych. Patrz Ropp, Shadid i Ober (2004) oraz Knoll, Chacon, Margolin i Mousseau (2003), aby zobaczyć perspektywę fizyków obliczeniowych na ten temat. Analiza matematyka pod względem warunków zamówienia znajduje się w książce Hairera i Wannera na temat sztywnej ODE (metody Rosenbrock-W są metodą implictu liniowo domniemaną), Kennedy i Carpenter (2003) dla nieliniowo domniemanego „dodatku” IMEX Runge-Kutta, oraz strona Emila Constantinescu, aby uzyskać najnowsze metody IMEX.

Ogólnie rzecz biorąc, metody IMEX mają więcej warunków zamówienia niż same ukryte metody ukryte i jawne. Pary metod IMEX można projektować z pożądaną stabilnością liniową i nieliniową, tak aby spełniały wszystkie warunki rzędu aż do kolejności projektowania metody. Spełnienie wszystkich warunków zamówienia pozwoli zachować asymptotyczny błąd podziału w tej samej skali, co błąd w każdym schemacie osobno. Nie mówi nic o reżimie przed asymptotycznym (duże kroki czasowe / wymaganie niskiej dokładności), ale rzadko jest bardziej rygorystyczne niż rozdzielczość każdej części osobno. W każdym razie błąd podziału jest widoczny dla wbudowanego estymatora błędów (w przypadku korzystania z adaptacyjnej kontroli błędów).

PETSc ma wiele metod IMEX z rodzin Rosenbrock-W i dodatków Runge-Kutta , i będzie miał ekstrapolację i liniowy wieloetapowy IMEX w naszym następnym wydaniu.

Oświadczenie: Napisałem wiele ze wsparcia integracji czasu PETSc i współpracuję z Emilem (link powyżej).

— Jed Brown
źródło

Z pewnością podchodzę do tego z perspektywy fizyki, więc wszystkie szczegóły techniczne zajmują mi trochę czasu, ponieważ nie znam wielu terminów. Jestem eksperymentalistą! Czy mógłbyś wyjaśnić nieco więcej na temat warunków zamówienia? IMEX czy te wieloetapowe metody są wymienione przez Jana?

— boyfarrell

Warunki porządku są zależnościami między współczynnikami metod ODE (np. Wpisy w tablicy Butchera dla metod Runge-Kutta), które muszą być spełnione, aby mieć porządek dokładności. Warunki zamówienia są omówione w dowolnej książce lub pracy, która projektuje metody integracji ODE, ale w zasadzie sprowadza się to do wielokrotnego stosowania pochodnych i dopasowywania warunków w rozwinięciu Taylora. Liczba warunków rzędu szybko rośnie dla metod wyższego rzędu, dlatego trudno jest zaprojektować metody wyższego rzędu. Bariery ustala się, pokazując, że warunki zamówienia są wzajemnie niezgodne.

— Jed Brown