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

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 ,

$\frac{\partial u}{\partial t} + \nabla \cdot \left( \boldsymbol{v} u - D\nabla u \right) = f$

gdzie $D$ to współczynnik dyfuzji substancji $u$ a $\boldsymbol{v}$ to prędkość.

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

$\frac{\partial u}{\partial t} = \underbrace{D\frac{\partial^2 u}{\partial x^2}}_{\textrm{Diffusion}} + \underbrace{\boldsymbol{v}\frac{\partial u}{\partial x}}_{\textrm{Advection (convection)}} + \underbrace{f(x,t)}_{\textrm{Reaction}}$

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

$\frac{u_{j}^{n+1} - u_{j}^{n}}{\Delta t} = D \left[ \frac{1 - \beta}{(\Delta x)^2} \left( u_{j-1}^{n} - 2u_{j}^{n} + u_{j+1}^{n} \right) + \frac{\beta}{(\Delta x)^2} \left( u_{j-1}^{n+1} - 2u_{j}^{n+1} + u_{j+1}^{n+1} \right) \right] + \boldsymbol{v} \left[ \frac{1-\alpha}{2\Delta x} \left( u_{j+1}^{n} - u_{j-1}^{n} \right) + \frac{\alpha}{2\Delta x} \left( u_{j+1}^{n+1} - u_{j-1}^{n+1} \right) \right] + f(x,t)$

Zwróć uwagę na warunki $\alpha$ i $\beta$ . Umożliwia to przełączanie schematu między:

• $\beta=\alpha=1/2$ Crank-Niscolson,
• $\beta=\alpha=1$ jest w pełni niejawny
• $\beta=\alpha=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ł?

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.$\alpha=\beta=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.$\alpha=\beta=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 .$\alpha\ne\beta$

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$$\alpha=\beta=1$$D$$\alpha=\beta=1/2$$D$

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ść.$(u_{j+1}-u_{j-1})/2\Delta t$$(u_j-u_{j-1})/\Delta 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 .$e^{i(k j\Delta x-\omega n\Delta t)}$$k$$\omega$$\xi^n e^{ikj\Delta x}$$|\xi|\le 1$

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

$\frac{u^{n+1}_j-u^n_j}{\Delta t} = D\frac{u^{n+1}_{j-1}-2u^{n+1}_j+u^{n+1}_{j+1}}{\Delta x^2}$

Jeśli podstawimy falę płaską, a następnie podzielimy przez i , otrzymamy równanie$\xi^n$$e^{ikj\Delta x}$

$\frac{\xi-1}{\Delta t} = D\frac{e^{-ik\Delta x}-2+e^{ik\Delta x}}{\Delta x^2}\xi$

Wyczyść to teraz trochę, a otrzymamy:

$\xi = \frac{1}{1+\frac{2D\Delta t}{\Delta x^2}\left(1-\cos k\Delta x\right)}$ .

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

$\frac{u_j^{n+1}-u_j^n}{\Delta t} = v\frac{u_{j-1}^n-u_{j+1}^n}{2\Delta 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.$\xi$$|\xi|^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.$u_j-u_{j-1}$$v > 0$$u_j-u_{j+1}$$v < 0$$V\Delta t/\Delta x \le 1$