Rozpatrzmy następujące zagadnienie brzegowe:

dla nieujemnej stałej c, ustalonego odcinka a,b i znanych wartości α,β.

Na podstawie tego modelowego zadania opiszemy ideę metody różnic skończonych (MRS).

Przyjmijmy następujące oznaczenia na różnicę skończoną w przód (ang. forward finite difference) i różnicę skończoną w tył (ang. backward finite difference) :

dla h>0. Będziemy często opuszczali dolny indeks h, jeśli h będzie ustalone.

Dla ustalonego kroku h>0 rozważmy następującą aproksymację drugiej pochodnej:

Nietrudno zauważyć, że jeśli funkcja u jest klasy C4 w otoczeniu x to:

co pozostawiamy jako zadanie, zob. ćwiczenie 7.1.

Wprowadzając siatkę (ang. mesh), czyli zbiór dyskretny dla h=b-a/N :

z Ωh=a+k*hk=1N-1 i ∂⁡Ωh=a+k*hk=0,N=a,b, możemy zdefiniować następujące zadanie dyskretne: znaleźć uh funkcję określoną na siatce Ω¯h taką, że

Przyjmujemy, że g:∂⁡Ω→R przyjmuje wartości g⁢a=α i g⁢b=β.

Możemy przypuszczać, że uh będzie aproksymowało w jakimś sensie u rozwiązanie zadania wyjściowego w punktach siatki.

Jak porównać uh określone na zbiorze dyskretnym z u określonym na całym domknięciu obszaru? Możemy porównać te funkcje w punktach siatki licząc błąd:

dla rh⁢u funkcji obcięcia (ang. restriction) określonej na siatce, przyjmującej wartości u w punktach siatki, ale możemy też badać dyskretną normę Euklidesową Lh2, czyli pierwiastek z sumy kwadratów błędów w punktach siatki przeskalowanych przez h:

dla xk=a+k*h. Popatrzmy na wykres rozwiązania zadania dyskretnego dla c=0,a=0,b=π i f⁢x=sin⁡x, dla którego rozwiązaniem jest u⁢x=sin⁡x dla sześciu punktów i dla dwudziestu punktów, por. rysunki 7.1 i 7.2. Widać, że dla dwudziestu punktów siatki rozwiązanie przybliżone pokrywa się z rozwiązaniem dokładnym w punktach siatki.

Będziemy badać błąd w normach ∥⋅∥∞,h i ∥⋅∥0,h dla połowionych kroków, tzn. dla 2-k⁢h0 dla ustalonego h0=π/10. Wyniki w tabeli 7.1 sugerują, że błędy dyskretne w obu normach są rzędu dwa, tzn. że rh⁢u-uh∞,h=O⁢h2 i rh⁢u-uh0,h=O⁢h2. W kolejnych rozdziałach wyjaśnimy dlaczego tak jest. Jak się okaże wyniki eksperymentu są zgodne z oszacowaniami otrzymanymi teoretycznie.

Przyjmując oznaczenie uk=uh⁢xk=uh⁢a+k*h otrzymujemy następujący układ równań liniowych:

Wstawiając u0=g⁢a i uN=g⁢b do układu równań otrzymujemy następujący układ równań:

czyli układ z macierzą trójdiagonalną, który można rozwiązać np. metodą przeganiania (wersja rozkładu LU dla macierzy trójdiagonalnej) kosztem O⁢N dla N≤106, czy nawet większych N (w zależności od dostępnego komputera).

W tym rozdziale rozpatrzymy modelowe zadanie eliptyczne dwuwymiarowe na kwadracie jednostkowym Ω¯=0,12. Zadanie polega na znalezieniu u∈C2⁢Ω∩C⁢Ω¯ takiego, że

gdzie △⁢u=∑k=12∂2⁡u∂⁡xk2, c jest ustaloną nieujemną stałą, f - to funkcja ciągła na Ω, a g - to funkcja ciągła na ∂⁡Ω. Zakładamy, że istnieje jednoznaczne rozwiązanie (7.8).

Dla ustalonego kroku h>0 rozważmy następującą aproksymację drugiej pochodnej cząstkowej (por. (7.3)):

gdzie e→k jest k-tym wersorem.

Wprowadzamy siatkę (zbiór dyskretny) dla h=1/N:

i jej odpowiednie podzbiory Ωh=k*h,l*hk,l=1N-1⊂Ω i ∂⁡Ωh=k*h,l*hk,l=0,N⊂∂⁡Ω.

Definiujemy następujące zadanie dyskretne: chcemy znaleźć uh funkcję określoną na siatce Ω¯h taką, że

Tak, jak w przypadku jednowymiarowym, możemy porównywać błąd w punktach siatki w dyskretnej normie maksimum:

dla rh⁢u zdefiniowanego analogicznie, tzn. funkcji określonej na siatce przyjmującej wartości u w punktach siatki lub w dyskretnej normie Lh2:

dla xk,l=k*h,l*h.

W Tabeli 7.2 podane są wyniki obliczeń dla dyskretyzacji (7.11) zadania (7.8) z c=0 ze znanym rozwiązaniem u⁢x=sin⁡x⁢cos⁡x i odpowiednio dobranymi f⁢x,y=2⁢sin⁡x⁢cos⁡x i g⁢x,y=u⁢x,y dla x na brzegu kwadratu. Stosunki norm dyskretnych dla danego h względem błędu dla 2*h sugerują, że w tym przypadku widzimy rząd zbieżności kwadratowej w obu normach, tzn. że rh⁢u-uh∞,h=O⁢h2 i rh⁢u-uh0,h=O⁢h2, tak samo jak w przypadku jednowymiarowym. Oczywiście obliczenia czyli rozwiązywanie odpowiedniego układu równań liniowych jest teraz bardziej kosztowne, jako że np. dla N=160 czyli h=1/160, otrzymujemy układ równań liniowych z M=1592=25281 niewiadomymi i z macierzą o M2=6.3912*108 elementach. Jednak macierz jest pasmowa - o paśmie szerokości N-1 i o około 4*M, czyli ma 105 niezerowych elementów. Zatem możemy jeszcze zastosować specjalne, bezpośrednie metody rozwiązywania układów równań liniowych, takie jak odpowiednia wersja rozkładu L⁢U, por. np. [9], czy - jeśli trzymamy tę macierz w formacie rzadkim, np. spakowanych kolumn, czy wierszy - to możemy zastosować jakąś bezpośrednią metodę dla macierzy rzadkich, np. metodę frontalną, por. [8]. Warto zauważyć, że w tym szczególnym przypadku obszaru Ω=0,12 znamy wartości i wektory własne tej macierzy i możemy zastosować specjalne metody rozwiązywania tego układu równań z wykorzystaniem algorytmu szybkiej transformaty Fouriera (FFT) (ang. Fast Fourier Transform), por. [10].

Moglibyśmy też rozważyć siatkę Ω¯h z różnymi rozmiarami hk k=1,2 względem obu osi tzn.:

dla h1=1/N i h2=1/M i jej odpowiednie podzbiory Ωh=Ω¯h∩Ω i ∂⁡Ωh=Ω¯h∩∂⁡Ω.

Wtedy oczywiście otrzymalibyśmy trochę inny układ równań, ale generalnie o tych samych właściwościach.