Czy Crank-Nicolson jest stabilnym schematem dyskretyzacji równania reakcji-dyfuzji-doradztwa (konwekcji)?

26

Nie znam zbyt dobrze wspólnych schematów dyskretyzacji dla PDE. Wiem, że Crank-Nicolson jest popularnym schematem dyskretyzacji równania dyfuzji. Czy jest to również dobry wybór na okres doradztwa?

Interesuje mnie rozwiązanie równania Reakcja-Dyfuzja-Doradztwo ,

ut+(vuDu)=f

gdzie D to współczynnik dyfuzji substancji u a v to prędkość.

Dla mojej konkretnej aplikacji równanie można zapisać w postaci,

ut=D2ux2Diffusion+vuxAdvection (convection)+f(x,t)Reaction

Oto schemat Crank-Nicolson, który zastosowałem,

ujn+1ujnΔt=D[1β(Δx)2(uj1n2ujn+uj+1n)+β(Δx)2(uj1n+12ujn+1+uj+1n+1)]+v[1α2Δx(uj+1nuj1n)+α2Δx(uj+1n+1uj1n+1)]+f(x,t)

Zwróć uwagę na warunki α i β . Umożliwia to przełączanie schematu między:

  • β=α=1/2 Crank-Niscolson,
  • β=α=1 jest w pełni niejawny
  • β=α=0 jest w pełni jawny

Wartości mogą być różne, co pozwala, aby wyrazem dyfuzji był Crank-Nicolson, a słowem wskazującym coś innego. Jakie jest najbardziej stabilne podejście, co byś polecił?

boyfarrell
źródło

Odpowiedzi:

15

To dobrze sformułowane pytanie i bardzo przydatna rzecz do zrozumienia. Korrok słusznie odsyła cię do analizy von Neumanna i książki LeVeque. Mogę dodać do tego trochę więcej. Chciałbym napisać szczegółową odpowiedź, ale w tej chwili mam czas tylko na krótką:

Z otrzymujesz metodę, która jest absolutnie stabilna dla dowolnie dużych rozmiarów kroków, a także dokładność drugiego rzędu. Jednak metoda ta nie jest stabilna na L , więc bardzo wysokie częstotliwości nie będą tłumione, co jest niefizyczne.α=β=1/2

Z otrzymujesz metodę, która jest również bezwarunkowo stabilna, ale dokładna tylko pierwszego rzędu. Ta metoda jest bardzo rozpraszająca. Jest stabilny na L.α=β=1

Jeśli weźmiesz , twoja metoda może być rozumiana jako zastosowanie addytywnej metody Runge-Kutty do półdyskretyzacji z różnicą centrowaną. Analiza stabilności i dokładności takich metod jest znacznie bardziej skomplikowana. Bardzo ładny artykuł na temat takich metod jest tutaj .αβ

To, które podejście polecić, zależy w dużej mierze od wielkości , rodzaju początkowych danych, z którymi masz do czynienia, oraz dokładności, której szukasz. Jeśli akceptowalna jest bardzo niska dokładność, wówczas jest bardzo solidnym podejściem. Jeśli jest umiarkowane lub duże, problem jest zdominowany przez dyfuzję i bardzo sztywny; zazwyczaj daje dobre wyniki. Jeśli jest bardzo małe, może być korzystne użycie jawnej metody i podniesienie wyższego rzędu dla terminów konwekcyjnych.Dα=β=1Dα=β=1/2D

David Ketcheson
źródło
Bardzo wnikliwe odpowiedzi, dziękuję! Czy istnieje sposób zdefiniowania różnych reżimów dominacji dyfuzji i dominacji poradnictwa? Innym niż porównanie wielkości terminów? Na przykład, porównując tylko współczynniki? Jakie jest znaczenie terminu „stabilność L”. Wszyscy polecają tę książkę, muszę ją kupić!
boyfarrell
Podane przeze mnie kryterium dotyczy tylko współczynników. W skrócie, stabilność L oznacza, że ​​wysokie częstotliwości będą silnie tłumione.
David Ketcheson
Więc kiedy jest funkcją płynną (jak w tym sensie, że nie ma komponentów Fouriera o wysokiej częstotliwości) Crank-Nicolson jest dobrym wyborem. Jeśli jednak ma ostre krawędzie, to jest dobrym wyborem. u(x)u(x)β=1
boyfarrell,
To rozsądne, choć bardzo szorstkie uogólnienie. Te wybory zadziałają przynajmniej, jeśli nie potrzebujesz dużej dokładności.
David Ketcheson
10

Ogólnie rzecz biorąc, będziesz chciał zastosować metodę niejawną dla równań parabolicznych (część dyfuzyjna) - wyraźne schematy parabolicznego PDE muszą mieć bardzo krótki czas, aby być stabilnym. I odwrotnie, dla części hiperbolicznej (porady) będziesz potrzebować jawnej metody, ponieważ jest ona tańsza i nie zakłóca symetrii układu liniowego, którą musisz rozwiązać za pomocą niejawnego schematu dyfuzji. W takim przypadku chcesz uniknąć wyśrodkowanych różnic, takich jak i przełącz się na różnice jednostronne ze względu na stabilność.(uj+1uj1)/2Δt(ujuj1)/Δt

Proponuję spojrzeć na książkę Randy'ego Leveque'a lub książkę Dale'a Durrana „Analiza stabilności von Neumanna”. Jest to ogólne podejście do ustalania stabilności programu dyskretyzacji, pod warunkiem, że masz okresowe warunki brzegowe. (Jest to także dobry artykuł wiki tutaj ).

Podstawową ideą jest założenie, że dyskretne przybliżenie można zapisać jako sumę fal płaskich , gdzie to liczba fal, a częstotliwość. Wpychasz falę samolotu w przybliżenie do PDE i modlisz się, żeby nie wybuchła. Możemy przepisać falę płaską jako i chcemy się upewnić, że .ei(kjΔxωnΔt)kωξneikjΔx|ξ|1

Jako przykład rozważmy równanie dyfuzji zwykłej z całkowicie niejawnym różnicowaniem:

ujn+1ujnΔt=Duj1n+12ujn+1+uj+1n+1Δx2

Jeśli podstawimy falę płaską, a następnie podzielimy przez i , otrzymamy równanieξneikjΔx

ξ1Δt=DeikΔx2+eikΔxΔx2ξ

Wyczyść to teraz trochę, a otrzymamy:

ξ=11+2DΔtΔx2(1coskΔx) .

Jest to zawsze mniej niż jeden, więc masz pewność. Spróbuj zastosować to do jawnego, wyśrodkowanego schematu równania porady:

ujn+1ujnΔt=vuj1nuj+1n2Δx

i zobaczyć, co otrzymujemy. (Tym razem będzie to część wyobrażona.) Przekonasz się, że , co jest smutnym czasem. Stąd moje upomnienie, że go nie używasz. Jeśli możesz to zrobić, nie powinieneś mieć większych problemów ze znalezieniem stabilnego schematu dla pełnego równania rada-dyfuzja.ξ|ξ|2>1

To powiedziawszy, użyłbym w pełni niejawnego schematu dla części dyfuzyjnej. Zmień różnicowanie w części na jeśli i jeśli i wybierz aby . (Jest to warunek Courant-Friedrichs-Lewy .) Jest dokładny tylko pierwszego rzędu, więc jeśli chcesz, możesz sprawdzić schematy dyskretyzacji wyższego rzędu.ujuj1v>0ujuj+1v<0VΔt/Δx1

Daniel Shapero
źródło
To bardzo szczegółowa odpowiedź, dziękuję.
boyfarrell
Ta odpowiedź uwzględnia tylko dyskretyzacje oparte na metodach Eulera do przodu i do tyłu w czasie. Pytanie dotyczy Cranka-Nicholsona.
David Ketcheson