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

Rozpatrzmy siatkę równomierną na półpłaszczyźnie 0,∞×R z krokiem przestrzennym h>0 i czasowym τ>0:

dla

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

lub

W tym rozdziale zakładamy, że spełniony jest warunek początkowy uk0=u0⁢xk dla k∈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

lub równoważnie biorąc λ=τh:

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

czyli

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

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

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 ukn w niestabilnym schemacie (15.5) zastępujemy średnią z uk+1n i uk-1n i otrzymujemy:

Kolejny schemat Laxa-Wendroffa wyprowadza się z rozwinięcia rozwiązania w momencie t=tn w szereg Taylora:

Korzystamy następnie z równania:

i kolejnego równania otrzymanego z poprzedniego:

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

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

i w końcu otrzymujemy schemat Laxa-Wendroffa:

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:

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:

gdzie Hk+1/2n=H⁢ukn,uk+1n 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:

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

przyjmując oznaczenie Fk+1/2n=H⁢F⁢ukn,F⁢ukn+1 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:

gdzie Λ=diag⁢λ1,…,λ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 Ukn=u→⁢tn,xk w punktach siatki, szukamy Wkn=C⁢Ukn, czyli stosujemy schematy do równoważnego równania (w→=C⁢u→):

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

Zatem możemy zastosować np. schemat upwind lub inny - niezależnie do każdej składowej - i otrzymać przybliżone rozwiązanie wj⁢tn,xk.

Znając wartości Wkn możemy wrócić do wyjściowych zmiennych i obliczyć Ukn rozwiązując układ równań

czyli przybliżone rozwiązanie u→⁢tn,xk.

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 λj i kolumny C. Następnie stosujemy wybrany schemat obliczając wartości Wkn 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 Ukn.

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 un=uknk∈Z i załóżmy, że dwupoziomowy (względem czasu) schemat różnicowy możemy opisać jako

lub w punkcie siatki xk

ze znanym warunkiem początkowym u0, tzn. uk0=u0⁢xk. Stabilność w pewnej normie ∥⋅∥h na odcinku czasu 0,T oznacza, że istnieją stałe τ0,C>0 takie, że dla czasów tn=n*τ<T jeśli 0<h,τ<τ0:

Często stosowaną normą jest norma będąca aproksymacją normy L1⁢R, czyli

Jeśli istnieją stałe τ0,β>0 takie, że dla czasów tn=n*τ<T jeśli 0<h,τ<τ0 zachodzi:

to

dla dowolnych n takich, że tn=n*τ<T, co oznacza stabilność schematu w normie ∥⋅∥h.

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

spełnia

dla u rozwiązania wyjściowego równania, dowolnego 0<h≤τ0 i t>0. Tutaj Aτ⁢u⁢t;k jest zdefiniowane dla tego schematu jak Aτ⁢un;k zastępując ukn przez u⁢t,xk.

Jeśli dla pewnej stałej C>0 i 0<h,τ≤τ0 zachodzi oszacowanie:

to powiemy, że rząd aproksymacji schematu wynosi q1 po czasie i q2 po przestrzeni. A jeśli zachodzi stała zależność τ od h np. liniowa τ=κ*h dla stałej κ, to mówimy, że rząd aproksymacji schematu wynosi q=min⁡q1,q2, 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:

dla h,τ→0 wtedy i tylko wtedy, jeśli jest stabilny i zgodny.

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.

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

gdzie γ∈C i α∈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 γ>1, to oczywiście schemat stabilny być nie może, a jeśli wszystkie takie rozwiązania dla dowolnych α muszą spełniać γ<1 - ewentualnie przy pewnych warunkach na τ i h - to schemat ma szanse być stabilnym, czy - inaczej: jest stabilnym w klasie rozwiązań tej postaci.

Podsumowując: wstawimy ukn=γn⁢ei⁢α⁢k do schematu i wyliczamy γ⁢α jeśli γ⁢α<1 dla dowolnego α, 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.