Zagadnienia

7.1 Modelowe zadanie jednowymiarowe
7.2 Modelowe zadanie dwuwymiarowe
- 7.2.1 Warunki brzegowe dla obszaru o skomplikowanej geometrii
7.3 Zadania

7. Metoda różnic skończonych dla równań eliptycznych drugiego rzędu

W tym rozdziale przedstawimy idee metody różnic skończonych na dwóch modelowych przykładach. Ze wszystkich metod przybliżonego rozwiązywania równań różniczkowych cząstkowych metoda ta wydaje się najbardziej intuicyjna w konstrukcji. W klasycznym sformułowaniu danego równania różniczkowego zamiast pochodnych rozpatrujemy ich przybliżenia. Na zadanej siatce rozpatrujemy przybliżenia pochodnych za pomocą różnic skończonych, czyli ilorazów różnicowych.

7.1. Modelowe zadanie jednowymiarowe

Rozpatrzmy następujące zagadnienie brzegowe:

$\displaystyle Lu(x)=-u^{{\prime\prime}}(x)+c*u(x)$	$\displaystyle=$	$\displaystyle f(x)\qquad x\in\Omega=(a,b),$	(7.1)
$\displaystyle u(a)$	$\displaystyle=$	$\displaystyle\alpha,$
$\displaystyle u(b)$	$\displaystyle=$	$\displaystyle\beta$

dla nieujemnej stałej $c$ , ustalonego odcinka $[a,b]$ i znanych wartości $\alpha,\beta$ .

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) :

	$\displaystyle\partial _{h}u(x)=h^{{-1}}(u(x+h)-u(x)),$				(7.2)
	$\displaystyle\overline{\partial}_{h}u(x)=h^{{-1}}(u(x)-u(x-h)).$

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:

$\partial\overline{\partial}_{h}u(x)=\partial\overline{\partial}u(x)=\overline{\partial}{\partial}u(x)=\frac{u(x-h)-2u(x)+u(x+h)}{h^{2}}.$

(7.3)

Nietrudno zauważyć, że jeśli funkcja $u$ jest klasy $C^{4}$ w otoczeniu $x$ to:

$u^{{\prime\prime}}(x)=\partial\overline{\partial}u(x)+O(h^{2}).$

(7.4)

co pozostawiamy jako zadanie, zob. ćwiczenie 7.1.

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

$\overline{\Omega}_{h}=\{ a+k*h\} _{{k=0}}^{N}$

z $\Omega _{h}=\{ a+k*h\} _{{k=1}}^{{N-1}}$ i $\partial\Omega _{h}=\{ a+k*h\} _{{k=0,N}}=\{ a,b\}$ , możemy zdefiniować następujące zadanie dyskretne: znaleźć $u_{h}$ funkcję określoną na siatce $\overline{\Omega}_{h}$ taką, że

	$\displaystyle-\partial\overline{\partial}u_{h}(x)+c*u_{h}(x)=f(x)\qquad x\in\Omega _{h}$				(7.5)
	$\displaystyle u(x)=g(x)\qquad x\in\partial\Omega _{h}$				(7.6)

Przyjmujemy, że $g:\partial\Omega\rightarrow\mathbb{R}$ przyjmuje wartości $g(a)=\alpha$ i $g(b)=\beta$ .

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

$\$

Rys. 7.1. Wykres rozwiązania dokładnego $u(x)=\sin(x)$ zadania (7.1) i rozwiązania przybliżonego (7.5) oznaczonego plusami dla siatki sześciopunktowej .

$\$

Rys. 7.2. Wykres rozwiązania dokładnego $u(x)=\sin(x)$ zadania (7.1) i rozwiązania przybliżonego (7.5) oznaczonego plusami dla dwudziestu punktów.

Jak porównać $u_{h}$ 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:

$\| u_{h}-r_{h}u\| _{{\infty,h}}=\max _{{x\in\overline{\Omega}_{h}}}|u_{h}(x)-u(x)|$

dla $r_{h}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ą $L^{2}_{h}$ , czyli pierwiastek z sumy kwadratów błędów w punktach siatki przeskalowanych przez $h$ :

$\| u_{h}-r_{h}u\| _{{0,h}}=\sqrt{h\sum _{{k=0}}^{N}|u_{h}(x_{k})-u(x_{k})|^{2}}$

dla $x_{k}=a+k*h$ . Popatrzmy na wykres rozwiązania zadania dyskretnego dla $c=0,a=0,b=\pi$ 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 $\|\cdot\| _{{\infty,h}}$ i $\|\cdot\| _{{0,h}}$ dla połowionych kroków, tzn. dla $2^{{-k}}h_{0}$ dla ustalonego $h_{0}=\pi/10$ . Wyniki w tabeli 7.1 sugerują, że błędy dyskretne w obu normach są rzędu dwa, tzn. że $\| r_{h}u-u_{h}\| _{{\infty,h}}=O(h^{2})$ i $\| r_{h}u-u_{h}\| _{{0,h}}=O(h^{2})$ . W kolejnych rozdziałach wyjaśnimy dlaczego tak jest. Jak się okaże wyniki eksperymentu są zgodne z oszacowaniami otrzymanymi teoretycznie.

$\begin{array}[]{|l|c|c|c|c|}\hline N&\| e_{h}\| _{{\infty,h}}&\| e_{h}\| _{{\infty,h}}/\| e_{{2h}}\| _{{\infty,2h}}&\| e_{h}\| _{{0,h}}&\| e_{h}\| _{{0,h}}/\| e_{{2h}}\| _{{0,2h}}\\ \hline 10&8.265e-03&&1.036e-02&\\ 20&2.059e-03&4.01e+00&2.580e-03&4.01e+00\\ 40&5.142e-04&4.00e+00&6.445e-04&4.00e+00\\ 80&1.285e-04&4.00e+00&1.611e-04&4.00e+00\\ 160&3.213e-05&4.00e+00&4.027e-05&4.00e+00\\ 320&8.032e-06&4.00e+00&1.007e-05&4.00e+00\\ 640&2.008e-06&4.00e+00&2.517e-06&4.00e+00\\ 1280&5.020e-07&4.00e+00&6.292e-07&4.00e+00\\ 2560&1.255e-07&4.00e+00&1.573e-07&4.00e+00\\ 5120&3.138e-08&4.00e+00&3.933e-08&4.00e+00\\ \hline\end{array}$

Tabela 7.1. Badanie błędu dyskretnego dla dyskretyzacji różnicami skończonymi zadania $-u^{{\prime\prime}}=\sin(x)$ dla $x\in[0,\pi]$ z $u(0)=u(\pi)=0$ dla którego znamy rozwiązanie $u(x)=\sin(x)$ . W kolumnie drugiej podajemy normę błędu $e_{h}=r_{h}u-u_{h}$ w dyskretnej normie maksimum, a w kolumnie czwartej w dyskretnej normie $L^{2}$ . W kolumnach trzeciej i piątej podajemy stosunek odpowiedniej normy błędu dla danego $h$ względem normy błędu dla $2*h$ .

Przyjmując oznaczenie $u_{k}=u_{h}(x_{k})=u_{h}(a+k*h)$ otrzymujemy następujący układ równań liniowych:

$\begin{array}[]{rcl}u_{0}&=&g(a),\\ \frac{1}{h^{2}}(-u_{{k-1}}+2*u_{k}-u_{{k-1}})+c*u_{k}&=&f(x_{k})=f_{k}\qquad k=1,\ldots,N-1,\\ u_{N}&=&g(b).\end{array}$

Wstawiając $u_{0}=g(a)$ i $u_{N}=g(b)$ do układu równań otrzymujemy następujący układ równań:

$\frac{1}{h^{2}}\left(\begin{array}[]{rrrrr}2+c*h^{2}&-1&0&\cdots&0\\ -1&2+c*h^{2}&-1&\vdots\\ 0&\ddots&\ddots&\ddots&\vdots\\ &&-1&2+c*h^{2}&-1\\ 0&\cdots&0&-1&2+c*h^{2}\end{array}\right)\left(\begin{array}[]{c}u_{1}\\ u_{2}\\ \vdots\\ u_{{N-2}}\\ u_{{N-1}}\end{array}\right)=\left(\begin{array}[]{c}f_{1}+\frac{1}{h^{2}}g(a)\\ u_{2}\\ \vdots\\ u_{{N-2}}\\ u_{{N-1}}+\frac{1}{h^{2}}g(b)\end{array}\right),$

(7.7)

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\leq 10^{6}$ , czy nawet większych $N$ (w zależności od dostępnego komputera).

7.2. Modelowe zadanie dwuwymiarowe

W tym rozdziale rozpatrzymy modelowe zadanie eliptyczne dwuwymiarowe na kwadracie jednostkowym $\overline{\Omega}=[0,1]^{2}$ . Zadanie polega na znalezieniu $u\in C^{2}(\Omega)\cap C(\overline{\Omega})$ takiego, że

	$\displaystyle Lu(x)=-\triangle u(x)+c*u(x)$	$\displaystyle=$	$\displaystyle f(x)\qquad x\in\Omega$		(7.8)
	$\displaystyle u(s)=g(s)\qquad s\in\partial\Omega,$

gdzie $\triangle u=\sum _{{k=1}}^{2}\frac{\partial^{2}u}{\partial x_{k}^{2}}$ , $c$ jest ustaloną nieujemną stałą, $f$ - to funkcja ciągła na $\Omega$ , a $g$ - to funkcja ciągła na $\partial\Omega$ . 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)):

$\partial\overline{\partial}_{{k,h}}u(x)=\partial\overline{\partial}_{k}u(x)=\overline{\partial}{\partial}_{k}u(x)=\frac{u(x-h\vec{e}_{k})-2u(x)+u(x+h\vec{e}_{k})}{h^{2}}.\qquad k=1,2,$

(7.9)

gdzie $\vec{e}_{k}$ jest $k$ -tym wersorem.

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

$\overline{\Omega}_{h}=\{(k*h,l*h)\} _{{k,l=0}}^{N}$

(7.10)

i jej odpowiednie podzbiory $\Omega _{h}=\{(k*h,l*h)\} _{{k,l=1}}^{{N-1}}\subset\Omega$ i $\partial\Omega _{h}=\{(k*h,l*h)\} _{{k,l=0,N}}\subset\partial\Omega$ .

Definiujemy następujące zadanie dyskretne: chcemy znaleźć $u_{h}$ funkcję określoną na siatce $\overline{\Omega}_{h}$ taką, że

$\left\{\begin{array}[]{rcl}-\sum _{{k=1,2}}\partial\overline{\partial}_{k}u_{h}(x)+c*u_{h}(x)&=&f(x)\qquad x\in\Omega _{h},\\ u(x)&=&g(x)\qquad x\in\partial\Omega _{h}\end{array}\right.$

(7.11)

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

$\| u_{h}-r_{h}u\| _{{\infty,h}}=\max _{{x\in\overline{\Omega}_{h}}}|u_{h}(x)-u(x)|$

dla $r_{h}u$ zdefiniowanego analogicznie, tzn. funkcji określonej na siatce przyjmującej wartości $u$ w punktach siatki lub w dyskretnej normie $L_{h}^{2}$ :

$\| u_{h}-r_{h}u\| _{{0,h}}=\sqrt{h^{2}\sum _{{k,l=0}}^{N}|u_{h}(x_{{k,l}})-u(x_{{k,l}})|^{2}}$

dla $x_{{k,l}}=(k*h,l*h)$ .

$\begin{array}[]{|l|c|c|c|c|}\hline N&\| z_{h}\| _{{\infty,h}}&\| z_{h}\| _{{\infty,h}}/\| z_{h}\| _{{\infty,2h}}&\| z_{h}\| _{{0,h}}&\| z_{h}\| _{{0,h}}/\| z_{h}\| _{{0,2h}}\\ \hline 10&5.211e-05&&2.789e-05&\\ 20&1.317e-05&3.96e+00&7.011e-06&3.98e+00\\ 40&3.298e-06&3.99e+00&1.755e-06&3.99e+00\\ 80&8.254e-07&4.00e+00&4.389e-07&4.00e+00\\ 160&2.064e-07&4.00e+00&1.097e-07&4.00e+00\\ \hline\end{array}$

Tabela 7.2. Badanie błędu dyskretnego dla dyskretyzacji różnicami skończonymi dwuwymiarowego zadania $-\triangle u=2*\sin(x)*\cos(x)$ na $x\in(0,1)^{2}$ , dla którego znamy rozwiązanie $u(x)=\sin(x)\cos(x)$ z odpowiednim warunkiem brzegowym Dirichleta na $\partial\Omega$ . W kolumnie drugiej podajemy normę błędu $z_{h}=r_{h}u-u_{h}$ w dyskretnej normie maksimum, a w kolumnie czwartej w dyskretnej normie $L^{2}$ . W kolumnach trzeciej i piątej podajemy stosunek odpowiednich norm dla danego $h$ względem normy błędu dla $2*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 $\| r_{h}u-u_{h}\| _{{\infty,h}}=O(h^{2})$ i $\| r_{h}u-u_{h}\| _{{0,h}}=O(h^{2})$ , 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=159^{2}=25281$ niewiadomymi i z macierzą o $M^{2}=6.3912*10^{8}$ elementach. Jednak macierz jest pasmowa - o paśmie szerokości $N-1$ i o około $4*M$ , czyli ma $10^{5}$ 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 $LU$ , 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 $\Omega=(0,1)^{2}$ 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ę $\overline{\Omega}_{h}$ z różnymi rozmiarami $h_{k}$ $k=1,2$ względem obu osi tzn.:

$\overline{\Omega}_{h}=\{(k*h_{1},l*h_{2})\} _{{k=0,\ldots,N;l=0,\ldots,M}}$

dla $h_{1}=1/N$ i $h_{2}=1/M$ i jej odpowiednie podzbiory $\Omega _{h}=\overline{\Omega}_{h}\cap\Omega$ i $\partial\Omega _{h}=\overline{\Omega}_{h}\cap\partial\Omega$ .

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

7.2.1. Warunki brzegowe dla obszaru o skomplikowanej geometrii

Zauważmy, że dla obszaru $\overline{\Omega}=[0,1]^{2}$ możemy tak dobrać kroki siatki, aby otrzymać w sposób naturalny punkty siatki leżące na $\partial\Omega$ . W przypadku obszarów o bardziej skomplikowanej geometrii, które nie są prostokątami, czy sumami prostokątów, często nie ma takiej możliwości. Taka sytuacja ma miejsce np. gdy $\Omega$ jest kołem.

$\$

Rys. 7.3. Siatka dla obszaru nieprostokątnego. Czarne punkty należą do $\partial\Omega _{h}$ .

$\$

Rys. 7.4. Siatka dla obszaru nieprostokątnego. Otoczenia siatkowe dla operatora $L_{h}=-\sum _{{k=1,2}}\partial\overline{\partial}_{k}$ .

Dla dowolnego $\Omega$ i naszej aproksymacji różnicowej Laplasjanu wprowadzamy pomocniczą siatkę na ${\mathbb{R}^{2}}$ postaci $S_{h}=\vec{x}_{0}+(k*h_{1},l*h_{2})_{{k,l\in{\mathbb{Z}}}}$ dla $x_{0}$ ustalonego punktu i $h_{1},h_{2}$ dodatnich kroków. Dla $x\in S_{h}$ i naszego operatora różnicowego $L_{h}u_{h}(x)=-\sum _{{k=1,2}}\partial\overline{\partial}_{k}u_{h}(x)+c*u_{h}(x)$ definiujemy otoczenie siatkowe (ang. mesh neighborhood) punktu $x$ jako podzbiór $S_{h}$ taki, aby $L_{h}u_{h}(x)$ było zdefiniowane poprzez wartości $u_{h}$ w $N_{h}(x)$ , czyli w tym przypadku otoczenie siatkowe $N_{h}(x)$ zawiera dany punkt $x$ i cztery sąsiednie punkty leżące na pionowej i poziomej prostej przechodzącej przez $x$ tzn. $N_{h}(x)=\{ x,x+h_{1}*\vec{e}_{1},x-h_{1}*\vec{e}_{1},x+h_{2}*\vec{e}_{2},x-h_{2}*\vec{e}_{2},\}$ , por. rysunek 7.4.

Wtedy definiujemy $\overline{\Omega}_{h}=S_{h}\cap\overline{\Omega}$ , a

$\Omega _{h}=\{ x\in S_{h}\cap\Omega:N_{h}(x)\subset\overline{\Omega}\}\subset\Omega\cap S_{h},$

czyli złożoną z takich punktów $\Omega$ , że ten punkt i sąsiednie punkty na osiach siatki należą do $\Omega$ . Wtedy $\partial\Omega _{h}=\overline{\Omega}_{h}\setminus\Omega _{h}$ , por. rysunek 7.3. (W przypadku bardziej skomplikowanych operatorów definicja $N_{h}(x)$ może być inna, a zatem i definicja $\Omega _{h}$ może być inna).

Następnie będziemy szukali funkcji zdefiniowanej na tej siatce tj. w $\overline{\Omega}_{h}=\Omega _{h}\cup\partial\Omega _{h}$ .

Zauważmy, że wtedy dla każdego $x\in\Omega _{h}$

$L_{h}u_{h}(x)=-\sum _{{k=1,2}}\partial\overline{\partial}_{k}u_{h}(x)+cu_{h}(x)=f_{h}(x)=f(x)$

jest poprawnie zdefiniowany. Pojawia się natomiast pytanie, jak postawić warunek brzegowy Dirichleta w punktach $\partial\Omega _{h}$ , które nie leżą na $\partial\Omega$ . Najprostszym rozwiązaniem dla $x\in\partial\Omega _{h}$ , który leży w $\Omega$ , jest postawienie w takim punkcie warunku:

$l_{h}u_{h}(x)=g(x_{s})=g_{h}(x),$

gdzie $x_{s}\in\partial\Omega$ taki, że odległość od $x$ jest nie większa niż $h$ . Dalej postępujemy jak w przypadku obszaru prostokątnego.

Rozwiązanie wyjściowego zadania $u$ jest określone we wszystkich punktach $\overline{\Omega}_{h}$ , ponieważ $\partial\Omega _{h}\subset\overline{\Omega}$ . Operator obcięcia definiujemy tak samo, tzn.:

$r_{h}u(x)=u(x)\qquad x\in\overline{\Omega}_{h}.$

Dla $u$ dostatecznie gładkiego otrzymujemy:

$\| L_{h}r_{h}u-f_{h}\| _{{\infty,h,\Omega _{h}}}:=\max _{{x\in\Omega _{h}}}|L_{h}r_{h}u(x)-f_{h}(x)|=O(h^{2})$

dla $h=\max\{ h_{1},h_{2}\}$ i tylko:

$\| l_{h}r_{h}u-g_{h}(x)\| _{{\infty,h,\partial\Omega _{h}}}:=\max _{{x\in\partial\Omega _{h}}}|l_{h}r_{h}u(x)-g_{h}(x)|=O(h).$

jeśli $\partial\Omega _{h}\not\subset\partial\Omega$ . Pokazanie tego pozostawiamy jako zadanie.

Oznacza to, że dokładność schematu na rozwiązaniu wyjściowego zagadnienia różniczkowego, czyli rząd aproksymacji schematu (który formalnie zostanie zdefiniowany w kolejnym rozdziale, por. definicja 8.4) wynosi jeden. Istnieją oczywiście metody podwyższania rzędu aproksymacji w tym przypadku poprzez odpowiednią interpolację warunków brzegowych z brzegu obszaru na punkty siatki $\partial\Omega _{h}$ leżące wewnątrz $\Omega$ .

7.3. Zadania

Ćwiczenie 7.1

Udowodnij, (7.4).

Ćwiczenie 7.2

Rozpatrzmy zadanie różniczkowe jednowymiarowe:

$-u^{{\prime\prime}}(x)+c(x)u(x)=f\qquad w\qquad[0,1]$

czyli równanie (7.5), ale z warunkiem brzegowym Neumanna : tzn. z $u^{{\prime}}(a)=g(a)$ i $u^{{\prime}}(b)=g(b)$ . Pokaż, że jeśli $c(x)\geq c_{0}>0$ to zadanie ma jednoznaczne rozwiązanie.

Rozważmy następującą dyskretyzację na siatce $\overline{\Omega}_{h}=\{ x_{k}\} _{{k=0,\ldots,N}}$ z $h=(b-a)/N$ i $x_{k}=k*h$ :

	$\displaystyle L_{h}u_{h}(x):=-\partial\overline{\partial}u_{h}(x)+c*u_{h}(x)=f(x)\qquad x\in\Omega _{h}$
	$\displaystyle l_{{h,a}}u_{h}:=\partial _{h}u(a)=g(a)\qquad l_{{h,b}}u_{h}=\overline{\partial}_{h}u(b)=g(b)$

Sprawdź, czy to zadanie dyskretne ma jednoznaczne rozwiązanie. Sformułuj je jako układ równań liniowych:

$A\vec{u}=\vec{F}$

dla $\vec{u}=(u_{0},\ldots,u_{N})^{T}$ , znajdując macierz $A$ i wektor prawej strony $\vec{F}$ . Czy ta macierz jest symetryczna, czy jest nieosobliwa? Czy jest trój-diagonalna? Jak rozwiązać powyższy układ możliwie małym kosztem?

Ćwiczenie 7.3

Rozpatrzmy zadanie i schemat z poprzedniego zadania. Zbadaj, dla jakiego możliwie dużego $p$ lokalny błąd schematu posiada rząd $p$ , czyli

$\max\{\max _{{x\in\Omega _{h}}}|L_{h}r_{h}u(x)-f(x)|,\max _{{s\in\{ a,b\}}}|l_{{h,s}}r_{h}u(s)-g(s)|\}=O(h^{p})$

Zakładamy dowolnie wysoką regularność rozwiązania zadania wyjściowego.

Ćwiczenie 7.4

Rozpatrzmy zadanie różniczkowe i schemat z ćwiczenia 7.2, ale z $c=\alpha=\beta=0$ . Pokaż, że zadanie nie ma jednoznacznego rozwiązania w ogólności, ale ma z dokładnością do stałej, o ile $\int _{a}^{b}f=0$ . Jaki warunek musi spełniać $f_{h}$ , aby zadanie dyskretne miało rozwiązanie? Sformułuj ten schemat jako układ równań liniowych, jak w ćwiczeniu 7.2. Pokaż, że jądro macierzy $A$ jest jednowymiarowe. Znajdź bazę jądra tej macierzy. Wykorzystując tę informację zaproponuj tanią (w sensie ilości operacji arytmetycznych) metodę znalezienia rozwiązania dyskretnego z dodatkowym warunkiem $u_{h}(a)=0$ .

Ćwiczenie 7.5

Analogicznie do przypadku jednowymiarowego, biorąc $x_{{kl}}=(k*h,l*h)$ ,

$u_{{k-1+(l-1)*(N-1)}}=u(x_{{kl}})$

i $\vec{u}=(u_{{k-1+(l-1)*(N-1)}})_{{k,l=1}}^{{N-1}}$ możemy zapisać zadanie dyskretne (7.11), jako układ równań liniowych $A_{h}\vec{u}=\vec{f}$ z macierzą $A_{h}$ i wektorem prawej strony $\vec{f}$ . Wyznacz tę macierz i ten wektor (por. (7.7) dla przypadku jednowymiarowego). Oblicz, ile elementów różnych od zera ma ta macierz. Policz, ile operacji arytmetycznych jest potrzebnych do rozwiązania tego układu równań liniowych z zastosowaniem metody Choleskiego, czyli rozkładu $A=L^{T}L$ dla $L$ macierzy dolnotrójkątnej w wersji dla macierzy pasmowych.

Ćwiczenie 7.6 (laboratoryjne)

Stwórz w octavie macierz z poprzedniego zadania dla $c=0$ , tzn. macierz układu równań powstałego z (7.11), jako macierz pełną i rzadką (można wykorzystać funkcję octave'a sparse()). Następnie rozwiąż ten układ dla wektora prawej strony odpowiadającego funkcji $f=-\triangle(u_{*})$ dla $u_{*}=x*(1-x)*y(1-y)$ .

Korzystając z narzędzi octave'a tic() i toc() sprawdź czas rozwiązywania dla różnych $N$ np. $N=20,80,320$ itp. dla obu typów macierzy. Czy różnica jest znacząca?
Zbadaj błąd dyskretny w normie maksimum (funkcja octave'a norm(wektor,'inf')) i w dyskretnej normie $L^{2}$ (np. używając funkcji octave'a norm(x,2)) między rozwiązaniem dokładnym $u_{*}$ , a rozwiązaniem dyskretnym dla $N=10,20,40,80,160$ . Czy eksperyment potwierdza teorię, tzn. czy dla podwojonych $N$ (czyli połowionych $h$ ) błąd maleje czterokrotnie?

Ćwiczenie 7.7 (laboratoryjne)

Rozwiąż równanie $-\frac{d^{2}u}{dx^{2}}+u=0$ z warunkami brzegowymi $u(0)=u(20)=1$ na $[0,20]$ przy pomocy metody różnic skończonych, tzn. rozwiąż zadanie (7.1) dla $a=0,b=20$ i $c(x)=1$ za pomocą schematu (7.5) dla $N=100,200$ . Policz normę maksimum w punktach siatki między dokładnymi wartościami rozwiązania $u^{*}(x)=(e^{{t-20}}+e^{{-t}})/(1+e^{{-20}})$ a rozwiązaniem dyskretnym $u_{h}$ zadania (7.5), Porównaj z wynikami metody strzałów zastosowanej do tego zadania tj. z ćwiczeniem 6.2.

Ćwiczenie 7.8 (częściowo laboratoryjne)

Rozpatrzmy przeskalowaną macierz z (7.7) (dla $c=0$ ):

$A_{{N-1}}=\left(\begin{array}[]{rrrrrr}2&-1&&&&\\ -1&2&-1&&&\\ &\ddots&\ddots&\ddots&&\\ &&-1&2&-1\\ &&&-1&2\\ \end{array}\right)\in\mathbb{R}^{{N-1,N-1}}$

Pokaż, że jej wartości własne to $\lambda _{k}=4*\sin^{2}(\frac{k*\pi}{2*N})$ dla $k=1,\ldots,N-1$ z odpowiednimi wektorami własnymi: $\vec{x}_{k}=\{\sin(\frac{k*\pi*j}{N})\} _{{n=1}}^{{N-1}}$ dla $k=1,\ldots,N-1$ . Oszacuj uwarunkowanie macierzy tzn. $\frac{\max _{k}\lambda _{k}}{\min _{k}\lambda _{k}}$ dla $N\gg 1$ .

Porównaj z wynikami otrzymanymi przy pomocy funkcji octave'a eig() i cond().

Wskazówka:

Korzystamy z wzorów trygonometrycznych: $\sin(a)+\sin(b)=2*\sin((a+b)/2)*\cos((a-b)/2)$ oraz $\cos(2*a)=1-2*\sin^{2}(a)$ .

Rozwiązanie:

Zauważmy, że biorąc $v_{{j,k}}=(v_{{k}})_{j}=\sin(\frac{k*\pi*j}{N})$ otrzymujemy, że $v_{{0,k}}=v_{{N,k}}=0$ . Zatem wystarczy sprawdzić, czy zachodzą równania $-v_{{j-1,k}}+2v_{{j,k}}-v_{{j+1,k}}=\lambda _{k}v_{{j,k}}$ dla $j=1,\ldots,N-1$ .

Biorąc $\alpha _{k}=k*\pi/N$ otrzymujemy:

$\displaystyle-v_{{j-1,k}}+2v_{{j,k}}-v_{{j+1,k}}$	$\displaystyle=$	$\displaystyle-(\sin((j-1)\alpha _{k})-\sin((j+1)\alpha _{k})+2\sin(j\alpha _{k})$
	$\displaystyle=$	$\displaystyle-2\sin(j\alpha _{k})\cos(\alpha _{k})+2\sin(j*\alpha _{k})$
	$\displaystyle=$	$\displaystyle 2(1-\cos(\alpha _{k}))\sin(j\alpha _{k})=4\sin^{2}(\alpha _{k}/2)\sin(j\alpha _{k})$
	$\displaystyle=$	$\displaystyle 4\sin^{2}(\alpha _{k}/2)v_{{j,k}}=4\sin^{2}(\frac{k\pi}{2N})v_{{j,k}}=\lambda _{k}*v_{{j,k}}.$