Zagadnienia

15.1 Schematy różnicowe dla równania skalarnego
15.2 Schematy dla równań nieliniowych lub układów równań
15.3 Stabilność, zgodność i zbieżność schematów
15.4 Metoda Fouriera badania stabilności
15.5 Zadania

15. Metody numeryczne rozwiązywania równań hiperbolicznych pierwszego rzędu

W tym rozdziale zajmiemy się metodami rozwiązywania równań hiperbolicznych pierwszego rzędu (por. rozdział 2.2.2). Przedstawimy konstrukcję kilku otwartych schematów różnicowych oraz podamy ideę zbieżności schematów za [26].

Konstrukcję schematów różnicowych przedstawimy dla modelowych równań, tzn. równań liniowych skalarnych, czyli będziemy szukali przybliżeń funkcji $u=u(t,x)$ takiej, że

$u_{t}+a(t,x)u_{x}=b(t,x),$

(15.1)

Zazwyczaj rozwiązania będą spełniały też warunek początkowy $u(0,x)=u_{0}(x)$ , przy czym najczęściej będziemy zakładać dla prostoty prezentacji, że $a$ jest stałą, a $b=0$ .

Pokażemy jak stosować te schematy dla równań nieliniowych i układów równań liniowych postaci:

$\vec{u}_{t}+A\vec{u}_{x}=0,$

(15.2)

gdzie $A$ - to stała macierz $m\times m$ diagonalizowalna w jakiejś bazie (ponieważ jest to układ hiperboliczny).

15.1. Schematy różnicowe dla równania skalarnego

W tym rozdziale zajmiemy się schematami różnicowymi dla równania skalarnego (15.1). Zakładamy, że $u$ spełnia dany warunek początkowy:

$u(0,x)=u_{0}(x)\qquad x\in\mathbb{R}.$

Rozpatrzmy siatkę równomierną na półpłaszczyźnie $[0,\infty)\times\mathbb{R}$ z krokiem przestrzennym $h>0$ i czasowym $\tau>0$ :

$(t_{n},x_{k})\qquad n\in\mathbb{N}\quad k\in\mathbb{Z}$

dla

$t_{n}=n*\tau,\qquad x_{k}=k*h.$

Możemy wtedy zdefiniować najprostszy otwarty schemat w sposób następujący:

$\tau^{{-1}}(u_{k}^{{n+1}}-u^{n}_{k})+h^{{-1}}a(u^{n}_{k}-u^{n}_{{k-1}})=0,$

lub

$\tau^{{-1}}(u_{k}^{{n+1}}-u^{n}_{k})+h^{{-1}}a(u^{n}_{{k+1}}-u^{n}_{k})=0.$

W tym rozdziale zakładamy, że spełniony jest warunek początkowy $u^{0}_{k}=u_{0}(x_{k})$ dla $k\in\mathbb{Z}$ .

Oba schematy są otwarte. Można postawić pytanie: który ma lepsze własności? Okazuje się, że stabilność tych schematów zależy od znaku parametru $a$ .

Schemat upwind definiujemy jako

$\mathrm{(Upwind)}\qquad\tau^{{-1}}(u_{k}^{{n+1}}-u^{k}_{n})=\left\{\begin{array}[]{lr}-h^{{-1}}a(u^{n}_{{k+1}}-u^{n}_{k})&\qquad a<0\\ -h^{{-1}}a(u^{n}_{k}-u^{n}_{{k-1}})&\qquad a>0\end{array}\right.$

(15.3)

lub równoważnie biorąc $\lambda=\frac{\tau}{h}$ :

$\mathrm{(Upwind)}\qquad(u_{k}^{{n+1}}-u^{k}_{n})+a\frac{\lambda}{2}(u^{n}_{{k+1}}-u^{n}_{{k-1}})=0.5\lambda|a|(u^{n}_{{k+1}}-2u^{n}_{k}+u^{n}_{{k-1}})$

(15.4)

Jeśli wprost dyskretyzujemy pochodną po przestrzeni za pomocą różnicy centralnej, to otrzymujemy następujący schemat:

$\tau^{{-1}}(u_{k}^{{n+1}}-u^{k}_{n})+a\frac{u^{n}_{{k+1}}-u^{n}_{{k+1}}}{2h}=0$

czyli

$u^{{n+1}}_{k}=u^{n}_{k}-a\frac{\lambda}{2}(u^{n}_{{k+1}}-u^{n}_{{k+1}}).$

(15.5)

Schemat ten niestety okazuje się być niestabilnym (wg definicji, która pojawi się później). Różni się on od poprzedniego schematu upwind (15.3) brakiem dodatkowego członu

$0.5\lambda|a|(u^{n}_{{k+1}}-2u^{n}_{k}+u^{n}_{{k-1}})=\tau*h*0.5|a|\frac{u^{n}_{{k+1}}-2u^{n}_{k}+u^{n}_{{k-1}}}{h^{2}},$

który aproksymuje, por. (7.3) i (7.4):

$\tau*h*0.5*|a|\frac{\partial^{2}u}{\partial x^{2}},$

Ten człon można traktować jako sztuczną numeryczną lepkość (ang. numerical dissipation or artificial viscosity) dodaną do niestabilnego schematu (15.5), dzięki której schemat upwind (15.3) jest stabilny.

Wyprowadzimy teraz kilka kolejnych otwartych schematów.

Rozważmy teraz schemat Laxa-Friedrichsa, w którym $u^{n}_{k}$ w niestabilnym schemacie (15.5) zastępujemy średnią z $u^{n}_{{k+1}}$ i $u^{n}_{{k-1}}$ i otrzymujemy:

$\mathrm{(Lax-Friedrichs)}\qquad u^{{n+1}}_{k}=0.5(u^{n}_{{k+1}}+u^{n}_{{k-1}})-a\frac{\lambda}{2}(u^{n}_{{k+1}}-u^{n}_{{k+1}}).$

(15.6)

Kolejny schemat Laxa-Wendroffa wyprowadza się z rozwinięcia rozwiązania w momencie $t=t_{n}$ w szereg Taylora:

$u(t_{n}+\tau,x_{k})=u(t_{n},x_{k})+\tau\frac{\partial u}{\partial t}(t_{n},x_{k})+0.5\tau^{2}\frac{\partial^{2}u}{\partial t^{2}}(t_{n},x_{k})+O(\tau^{3}).$

Korzystamy następnie z równania:

$\frac{\partial u}{\partial t}=-a\frac{\partial u}{\partial x}$

i kolejnego równania otrzymanego z poprzedniego:

$\frac{\partial^{2}u}{\partial t^{2}}=a^{2}\frac{\partial^{2}u}{\partial x^{2}}.$

W tych równaniach zastępujemy pochodną po przestrzeni ilorazem różnicowym centralnym, por. (4.2):

$\frac{\partial u}{\partial x}(t,x)\approx 0.5h^{{-1}}(u(t,x+h)-u(t,x-h)),$

a drugą pochodną po przestrzeni jej przybliżeniem różnicowym na trzech punktach, por. (7.3) i (7.4):

$\frac{\partial^{2}u}{\partial x^{2}}(t,x)\approx 0.5h^{{-2}}(u(t,x-h)-2*u(t,x)+u(t,x+h))$

i w końcu otrzymujemy schemat Laxa-Wendroffa:

$\mathrm{(Lax-Wendroff)}\qquad u^{{n+1}}_{k}=u^{n}_{k}-a\, 0.5\lambda(u^{n}_{{k+1}}-u^{n}_{{k-1}})+0.5a^{2}\lambda^{2}(u^{n}_{{k+1}}-2u^{n}_{k}+u^{n}_{{k-1}}).$

(15.7)

Można też rozważać schematy wielopoziomowe ze względu na czas, np. schemat skoku żaby (leap-frog), w którym pochodną po czasie dyskretyzujemy przez pochodną centralną tak samo, jak pochodną po przestrzeni. Otrzymujemy wówczas schemat trzypoziomowy:

$\mathrm{(Leap-frog)}\qquad u^{{n+1}}_{k}=u^{{n-1}}_{k}-a\lambda(u^{n}_{{k+1}}-u^{n}_{{k+1}}).$

(15.8)

15.2. Schematy dla równań nieliniowych lub układów równań

Zauważmy, że wszystkie dotąd rozważane dwupoziomowe schematy dla równań skalarnych, tzn. (15.7), (15.3), (15.5), można zapisać w zunifikowany sposób jako:

$u^{{n+1}}_{k}=u^{n}_{k}-\lambda(H^{n}_{{k+1/2}}-H^{n}_{{k-1/2}})$

gdzie $H^{n}_{{k+1/2}}=H(u^{n}_{k},u^{n}_{{k+1}})$ jest określana jako numeryczny strumień (ang. numerical flux).

W ten sposób możemy schematy dla równania (15.1) łatwo przenieść na przypadek nieliniowych równań hiperbolicznych postaci:

$u_{t}+\frac{\partial F(u)}{\partial x}=0$

gdzie $F$ - to dana funkcja. $F(u)$ nazywamy strumieniem dla funkcji $u$ . Wtedy każdy schemat można zapisać jako

$u^{{n+1}}_{k}=u^{n}_{k}-\lambda(F^{n}_{{k+1/2}}-F^{n}_{{k-1/2}})$

przyjmując oznaczenie $F^{n}_{{k+1/2}}=H(F(u^{n}_{k}),F(u^{{n+1}}_{k}))$ z $H$ numerycznym strumieniem wziętym z wyjściowego schematu. Zauważmy, że dla równania (15.1) zachodzi $F(u)=a*u$ .

Również schematy te można zastosować do układu (15.2). Wtedy widzimy, że dla nieosobliwej macierzy $C$ :

$A=C\Lambda C^{T},$

gdzie $\Lambda=\mathrm{diag}(\lambda _{1},\ldots,\lambda _{m})$ - to macierz diagonalna z wartościami własnymi $A$ na diagonali (ponieważ jest to układ równań hiperbolicznych). Wtedy możemy zamienić zmienne: zamiast szukać wartości rozwiązania $U^{n}_{k}=\vec{u}(t_{n},x_{k})$ w punktach siatki, szukamy $W^{n}_{k}=CU^{n}_{k}$ , czyli stosujemy schematy do równoważnego równania ( $\vec{w}=C\vec{u}$ ):

$\vec{w}_{t}+\Lambda\vec{w}_{x}=0.$

Proszę zauważyć, że to równanie jest układem $m$ niezależnych równań hiperbolicznych (15.1), tzn.:

$(w_{j})_{t}+\lambda _{j}(w_{j})_{x}=0\qquad j=1,\ldots,m.$

Zatem możemy zastosować np. schemat upwind lub inny - niezależnie do każdej składowej - i otrzymać przybliżone rozwiązanie $w_{j}(t_{n},x_{k})$ .

Znając wartości $W^{n}_{k}$ możemy wrócić do wyjściowych zmiennych i obliczyć $U^{n}_{k}$ rozwiązując układ równań

$CU^{n}_{k}=W^{n}_{k},$

czyli przybliżone rozwiązanie $\vec{u}(t_{n},x_{k})$ .

W praktyce - w pierwszym kroku przeprowadzamy obliczenia wstępne rozwiązując numerycznie zadanie własne dla macierzy $A$ , tzn. obliczając wartości i wektory własne $A$ , czyli $\lambda _{j}$ i kolumny $C$ . Następnie stosujemy wybrany schemat obliczając wartości $W^{n}_{k}$ dla punktów siatki (w praktyce musimy ograniczyć zakres $k$ i $n$ ). Na koniec rozwiązujemy układ równań z macierzą $C$ dla odpowiednich $k$ i $n$ otrzymując $U^{n}_{k}$ .

15.3. Stabilność, zgodność i zbieżność schematów

Do schematów różnicowych dla równań hiperbolicznych stosuje się ogólną teorię zbieżności Laxa-Richmyera schematów różnicowych, analogiczną do teorii zbieżności z rozdziału 8.1, por. rozdział 14.2 w [26]. Aby uzyskać oszacowanie błędu w pewnej normie dyskretnej, należy wykazać odpowiedni rząd aproksymacji schematu i stabilność w tej normie.

Przy przyjętych powyżej oznaczeniach przyjmijmy, że $\mathbf{u}^{n}=\{ u^{n}_{k}\} _{{k\in\mathbb{Z}}}$ i załóżmy, że dwupoziomowy (względem czasu) schemat różnicowy możemy opisać jako

$\mathbf{u}^{n}=A_{\tau}(\mathbf{u}^{{n-1}})\qquad n>0,$

(15.9)

lub w punkcie siatki $x_{k}$

$u^{n}_{k}=A_{\tau}(u^{{n-1}};k)\qquad n>0.$

ze znanym warunkiem początkowym $\mathbf{u}^{0}$ , tzn. $u^{0}_{k}=u_{0}(x_{k})$ . Stabilność w pewnej normie $\|\cdot\| _{h}$ na odcinku czasu $[0,T]$ oznacza, że istnieją stałe $\tau _{0},C>0$ takie, że dla czasów $t_{n}=n*\tau<T$ jeśli $0<h,\tau<\tau _{0}$ :

$\|\mathbf{u}^{n}\| _{h}\leq C_{T}\|\mathbf{u}^{0}\| _{h}\qquad n>0.$

Często stosowaną normą jest norma będąca aproksymacją normy $L^{1}(\mathbb{R})$ , czyli

$\|\mathbf{u}^{n}\| _{h}=\sum _{{k\in\mathbb{Z}}}h|u^{n}_{k}|.$

Jeśli istnieją stałe $\tau _{0},\beta>0$ takie, że dla czasów $t_{n}=n*\tau<T$ jeśli $0<h,\tau<\tau _{0}$ zachodzi:

$\| A_{\tau}\mathbf{u}\| _{h}\leq(1+\beta\tau)\|\mathbf{u}\| _{h}\qquad\forall\mathbf{u},$

$\|\mathbf{u}^{n}\| _{h}\leq(1+\beta\tau)^{n}\|\mathbf{u}^{0}\| _{h}\leq\exp(\beta T)\|\mathbf{u}^{0}\| _{h}\qquad n>0$

dla dowolnych $n$ takich, że $t_{n}=n*\tau<T$ , co oznacza stabilność schematu w normie $\|\cdot\| _{h}$ .

Z kolei zgodność schematu (aproksymacja schematem wyjściowego zadania; ang. consistency) oznacza, że schemat dyskretny aproksymuje wyjściowe równanie, tzn.

$(E_{\tau}(t))_{k}:=\tau^{{-1}}|u(t+\tau,x_{k})-A_{\tau}(\mathbf{u}(t);k)|$

spełnia

$\lim _{{\tau\rightarrow 0}}\| E_{\tau}(t)\| _{h}=0$

dla $u$ rozwiązania wyjściowego równania, dowolnego $0<h\leq\tau _{0}$ i $t>0$ . Tutaj $A_{\tau}(\mathbf{u}(t);k)$ jest zdefiniowane dla tego schematu jak $A_{\tau}(\mathbf{u}^{n};k)$ zastępując $u^{n}_{k}$ przez $u(t,x_{k})$ .

Jeśli dla pewnej stałej $C>0$ i $0<h,\tau\leq\tau _{0}$ zachodzi oszacowanie:

$\| E_{\tau}(t)\| _{h}\leq C\,[\tau^{{q_{1}}}+h^{{q_{2}}}]$

to powiemy, że rząd aproksymacji schematu wynosi $q_{1}$ po czasie i $q_{2}$ po przestrzeni. A jeśli zachodzi stała zależność $\tau$ od $h$ np. liniowa $\tau=\kappa*h$ dla stałej $\kappa$ , to mówimy, że rząd aproksymacji schematu wynosi $q=\min\{ q_{1},q_{2}\}$ , co jest zgodne z definicją z rozdziału 8.1.

Jak już wiemy, teoria Laxa mówi, że stabilność i aproksymacja (zgodność) dają zbieżność. Tak jest też w tym przypadku, tzn. teoria Laxa-Richmyera (por. [27]) mówi, że schemat jest zbieżny:

$\max _{{0\leq n\leq T/\tau}}\| u(t_{n},\cdot)-\mathbf{u}^{n}\| _{h}\rightarrow 0$

dla $h,\tau\rightarrow 0$ wtedy i tylko wtedy, jeśli jest stabilny i zgodny.

Przykład 15.1

Rozpatrzmy najprostszy schemat Laxa-Friedrichsa (15.6). Jeśli założymy, że

$|a\lambda|\leq 1$

(15.10)

to otrzymujemy:

	$\displaystyle\\|\mathbf{u}^{{n+1}}\\| _{h}$	$\displaystyle=$	$\displaystyle h\sum _{k}\|u^{{n+1}}_{j}\|\leq\frac{h}{2}\left[(1-\lambda a)\sum _{j}\|u^{n}_{{j+1}}\|+(1+\lambda a)\sum _{j}\|u^{n}_{{j-1}}\|\right]$
		$\displaystyle=$	$\displaystyle\frac{1}{2}\left[(1-\lambda a)\\|\mathbf{u}^{n}\\| _{h}+(1+\lambda a)\\|\mathbf{u}^{n}\\| _{h}\right]=\\|\mathbf{u}^{n}\\| _{h}.$

czyli, że schemat jest stabilny warunkowo przy założeniu $|a\lambda|\leq 1$ .

Analogicznie można pokazać, że przy założeniu, że spełniony jest warunek (15.10), schemat Laxa-Wendroffa (15.7) i schemat upwind (15.3) są stabilne. Widzimy, że schemat leap-frog (15.8) jest stabilny (w sensie analogicznej definicji odpowiedniej dla schematów trzypoziomowych) przy założeniu, że w warunku (15.10) zachodzi ostra nierówność.

Natomiast schemat (15.5) przy założeniach typu (15.10) nie jest stabilny.

Wykazanie, że wszystkie wymienione schematy są zgodne i zbadanie jakie mają rzędy pozostawiamy jako zadanie.

15.4. Metoda Fouriera badania stabilności

Metoda opisana w tym rozdziale służy badaniu stabilności schematów, choć - de facto - może tylko sprawdzić stabilność w sensie negatywnym. Tzn. metoda pozwala wykazać, że jakiś schemat nie jest stabilny.

Zakładamy, że rozpatrujemy schemat dwupoziomowy lub trzypoziomowy, np. dający się zapisać jako (15.9), i że szukamy jego rozwiązań w postaci

$u^{n}_{k}=\gamma^{n}\exp(i\alpha k),$

gdzie $\gamma\in\mathbb{C}$ i $\alpha\in\mathbb{R}$ są stałymi.

Metoda polega na wyznaczeniu warunków na te stałe w zależności od konkretnej postaci schematu. Jeśli się okaże, że istnieje rozwiązanie tej postaci z $|\gamma|>1$ , to oczywiście schemat stabilny być nie może, a jeśli wszystkie takie rozwiązania dla dowolnych $\alpha$ muszą spełniać $|\gamma|<1$ - ewentualnie przy pewnych warunkach na $\tau$ i $h$ - to schemat ma szanse być stabilnym, czy - inaczej: jest stabilnym w klasie rozwiązań tej postaci.

Podsumowując: wstawimy $u^{n}_{k}=\gamma^{n}e^{{i\alpha k}}$ do schematu i wyliczamy $\gamma(\alpha)$ jeśli $|\gamma(\alpha)|<1$ dla dowolnego $\alpha$ , to schemat uważamy za stabilny w sensie opisanym powyżej.

Zbadanie stabilności schematów (15.7), (15.3), (15.8), (15.5) przy pomocy tej metody pozostawiamy jako zadanie.

15.5. Zadania

Ćwiczenie 15.1

Zbadaj przy pomocy metody Fouriera stabilność schematów:

Laxa Wendroffa (15.7),
schematu upwind (15.3),
schematu Leap-frog (15.8),
schematu opartego na różnicy centralnej (15.5)

Ćwiczenie 15.2 (Laboratoryjne)

Zaimplementuj w octave schemat Laxa-Wendroffa dla równania $au_{{x}}=u_{{t}}$ dla $a=1,-1,100,-100$ , przyjmując, że znamy rozwiązanie początkowe na $[-1,1]$ i warunki brzegowe $u(t,-1)=u(t,1)=0$ . Zbadaj rząd metodą połowionych kroków, czyli dla $h$ i $h/2$ policz błędy w ustalonym punkcie i ich stosunek.

Treść automatycznie generowana z plików źródłowych LaTeXa za pomocą oprogramowania wykorzystującego LaTeXML.

	$\displaystyle\\|\mathbf{u}^{{n+1}}\\| _{h}$	$\displaystyle=$	$\displaystyle h\sum _{k}\|u^{{n+1}}_{j}\|\leq\frac{h}{2}\left[(1-\lambda a)\sum _{j}\|u^{n}_{{j+1}}\|+(1+\lambda a)\sum _{j}\|u^{n}_{{j-1}}\|\right]$
		$\displaystyle=$	$\displaystyle\frac{1}{2}\left[(1-\lambda a)\\|\mathbf{u}^{n}\\| _{h}+(1+\lambda a)\\|\mathbf{u}^{n}\\| _{h}\right]=\\|\mathbf{u}^{n}\\| _{h}.$

Numeryczne rozwiązywanie równań różniczkowych