Zagadnienia

8.1 Laplasjan: dyskretyzacja metodą różnic skończonych
8.2 Równanie dyfuzji
8.3 Uwagi i uzupełnienia

8. Stacjonarne równanie dyfuzji

W bardzo wielu praktycznych modelach (biologicznych, chemicznych, fizycznych i nawet ekonomicznych) należy rozwiązać równanie różniczkowe postaci

		$\displaystyle-\operatorname{div}(a(x)\nabla u(x))=f(x),\qquad x\in\Omega,$		(8.1)
		$\displaystyle u(x)=g(x)\qquad x\in\partial\Omega.$		(8.2)

Niewiadomą jest funkcja $u:\bar{\Omega}\rightarrow\mathbb{R}$ . W klasycznym modelu fizycznym dyfuzji, $u$ odpowiadałoby stężeniu jakiejś substancji, a współczynnik $a$ nazywalibyśmy współczynnikiem dyfuzji. Zadana funkcja $f$ byłaby interpretowana jako zewnętrzne źródło substancji, a warunek brzegowy typu Dirichleta oznaczałby, że na brzegu obszaru utrzymujemy zadane stężenie równe $g$ . Warto pamiętać, że równaniem tego samego lub podobnego typu można opisywać wiele innych sytuacji, na przykład rozchodzenie się ciepła w pewnym materiale.

Ze względu na powszechność występowania takich i podobnych modeli w zastosowaniach matematyki, trudno wyobrazić sobie, by absolwent matematyki stosowanej nie potrafił skutecznie zmierzyć się z takim zadaniem, choćby w najprostszej postaci — wykorzystując w tym celu choćby najprostsze metody. Dlatego w niniejszym rozdziale zajmiemy się szczegółowo rozwiązaniem takiego zadania w Octave.

W ogólnym przypadku, moglibyśmy mieć do czynienia z macierzą współczynników $a(x)$ , ale tutaj dla uproszczenia (zadanie i tak będzie obliczeniowo wystarczająco skomplikowane) przyjmiemy, że

Współczynnik dyfuzji jest dodatnią funkcją skalarną, $a:\bar{\Omega}\rightarrow\mathbb{R}_{+}$ i ograniczoną: istnieją stałe absolutne $m,M$ takie, że $0<m\leq|a(\cdot)|\leq M$ .
$\Omega$ jest kostką w $\mathbb{R}^{d}$ , przy czym $d\in\{ 1,2,3\}$ — czyli rozpatrujemy dyfuzję w co najwyżej trójwymiarowym obszarze o bardzo prostej geometrii. Co robić w przypadku, gdy $\Omega$ nie jest kostką, powiemy parę słów pod koniec rozdziału.

Zaczniemy więc od samodzielnej dyskretyzacji równania (8.1) w najprostszym przypadku.

8.1. Laplasjan: dyskretyzacja metodą różnic skończonych

Na początek przyjmijmy, że warunek brzegowy jest jednorodny: $g=0$ oraz $a(x)\equiv 1$ . Wówczas zadanie (8.1) upraszcza się do równania Poissona:

		$\displaystyle-\Delta u(x)=f(x)\qquad x\in\Omega,$		(8.3)
		$\displaystyle u(x)=0\qquad x\in\partial\Omega.$		(8.4)

Najprostszą, choć zazwyczaj najmniej skuteczną (przez swój brak wyrafinowania i dostosowania do faktycznego przebiegu rozwiązania) metodą dyskretyzacji takiego równania jest użycie równomiernej siatki na $\Omega$ i metody różnic skończonych.

8.1.1. Dyskretyzacja jednowymiarowego laplasjanu

Gdy $d=1$ , to $\Omega$ jest odcinkiem¹⁶Bez zmniejszenia ogólności założymy, że początek odcinka znajduje się w punkcie $x=0$ ., $\Omega=(0,X)$ , a (8.3) staje się równaniem różniczkowym zwyczajnym z warunkiem brzegowym:

		$\displaystyle-\dfrac{d^{2}}{dx^{2}}u(x)=f(x)\qquad\forall x\in(0,X),$		(8.5)
		$\displaystyle u(0)=u(X)=0.$		(8.6)

Aby znaleźć jego przybliżone rozwiązanie, możemy na przykład zastąpić równanie różniczkowe odpowiadającym mu równaniem różnicowym,

		$\displaystyle-\dfrac{u_{{i-1}}-2u_{i}+u_{{i+1}}}{h_{x}^{2}}=f_{i}\qquad\forall i=1,\ldots,N_{x},$		(8.7)
		$\displaystyle u_{0}=u_{{N_{x}+1}}=0,$		(8.8)

gdzie $u_{i}$ ma odpowiadać wartości rozwiązania w węźle dyskretyzacji $x_{i}$ , $u_{i}\approx u(x_{i})$ , oraz $x_{i}=i\cdot h_{x}$ , przy czym $h_{x}=X/(N_{x}+1)$ . Odpowiedzią na pytanie, jak dobrze (i czy w ogóle) takie rozwiązanie dyskretne aproksymuje rozwiązanie dokładne, zajmujemy się na wykładzie z Numerycznych równań różniczkowych; tutaj interesuje nas przede wszystkim aspekt praktyczny, czyli metoda uzyskania konkretnego przybliżenia. Najpierw zobaczmy więc, do jakiego zagadnienia algebraicznego prowadzi równanie różnicowe (8.7).

Jest to oczywiście równanie liniowe na współczynniki $U_{{N_{x}}}=(u_{1},\ldots,u_{{N_{x}}})^{T}$ ,

$T_{{N_{x}}}\, U_{{N_{x}}}=F_{{N_{x}}},$

gdzie

$T_{{N_{x}}}=\dfrac{1}{h_{x}^{2}}\begin{pmatrix}2&-1&&&\\ -1&2&-1&&\\ &\ddots&\ddots&\ddots&\\ &&-1&2&-1\\ &&&-1&2\end{pmatrix}_{{N_{x}\times N_{x}}}.$

(8.9)

Ponieważ interesują nas dobre przybliżenia (to znaczy przypadek, gdy $h_{x}$ jest małe — lub bardzo małe), znaczy to, że $N_{x}$ będzie (bardzo) duże. Nasza macierz jest macierzą rozrzedzoną, bo ma nieco mniej niż $3N_{x}$ niezerowych elementów. Jej współczynnik wypełnienia (density), czyli stosunek liczby niezerowych do liczby wszystkich elementów macierzy, rzędu $1/N_{x}$ , maleje ze wzrostem $N_{x}$ . Nie ma więc sensu pamiętanie wszystkich $N_{x}^{2}$ jej elementów — raczej, pamiętalibyśmy tylko jej niezerowe elementy. Do tego świetnie nadaje się format macierzy rozrzedzonych Octave. Dowolną macierz w takim formacie można utworzyć poleceniem sparse, ale ze względu na to, że nasza konkretna macierz ma wyjątkowo prostą strukturę, użyjemy polecenia

hx = X/(Nx+1);
Tx = spdiags( repmat([-1 2 -1]/(hx^2), Nx, 1), [-1,0,1], Nx, Nx);

gdzie wcześniej musimy zdefiniować X $=X$ Nx $=N_{x}$ . Komenda repmat klonuje daną macierz Nx razy.

Dzięki reprezentacji naszej macierzy w formacie rzadkim, macierz rozmiaru $N_{x}$ zajmie nam w pamięci tylko około $8\cdot 3\cdot N_{x}=24\, N_{x}$ bajtów (w Octave wszystkie liczby rzeczywiste są domyślnie reprezentowane w podwójnej precyzji, a więc na 64 bitach, czyli 8 bajtach) — a nie standardowe $8\, N_{x}^{2}$ bajtów w przypadku, gdybyśmy reprezentowali ją jako macierz pełną.

Ćwiczenie 8.1

Spróbuj zaimplementować obliczanie $T_{x}$ , zastępując w kodzie powyżej

[-1 2 -1]/(hx^2)

przez pozornie ,,równie dobrze wyglądające”

[-1 2 -1]/hx*hx

Czy dostajesz tę samą macierz Tx? To jeszcze jeden przyczynek do tego, ja ważne jest, by na różne sposoby zweryfikować poprawność konkretnej implementacji.

8.1.2. Rozwiązanie równania Poissona w obszarze jednowymiarowym

Otrzymana powyżej macierz jest trójdiagonalna (a przy tym: symetryczna i dodatnio określona), zatem właściwie zrealizowany algorytm eliminacji Gaussa (lub lepiej: rozkładu $LDL^{T}$ ) pozwoli nam wyznaczyć rozwiązanie układu $T_{{N_{x}}}U_{{N_{x}}}=F_{{N_{x}}}$ kosztem $O(N_{x})$ flopów — a więc, z dokładnością do stałej, optymalnym. Standardowy operator rozwiązywania równań w Octave — ,,\” — sprawdza na wstępie, czy macierz układu jest właśnie rozrzedzona, a jeśli tak — stosuje odpowiedni solver oparty na eliminacji Gaussa, w naszym przypadku w pełni uwzględniający trójdiagonalną strukturę. Dlatego rozwiązanie naszego układu skutecznie i szybko obliczymy komendą

Ux = Tx \ Fx;

Ćwiczenie 8.2

Wyznacz przybliżone rozwiązanie równania

		$\displaystyle-k\,\dfrac{d^{2}}{dx^{2}}u(x)=f(x),\qquad x\in(0,2),$
		$\displaystyle u(0)=u(2)=0,$

gdzie $k=2.45$ oraz

$f(x)=\begin{cases}1/(2-x),\qquad x\in(0,1],\\ 1/x^{5},\qquad x\in[1,2).\end{cases}$

Rozwiązanie: [pokaż]

Dzieląc równanie przez stałą $k$ dostajemy (8.20) dla prawej strony równej $f(x)/k$ . Ponieważ $f(x)$ załamuje się w $x=1$ , weźmiemy siatkę, w której jeden z węzłów wypadnie akurat w $x=1$ . W przypadku siatki jednorodnej zagwarantujemy to sobie, biorąc nieparzyste $N_{x}$ . Wstępnie wydaje się, że biorąc $N_{x}=99\, X+1=199$ , uzyskamy dostatecznie dużą rozdzielczość schematu (wszak wówczas $h_{x}=10^{{-2}}$ ). Potem postaramy się sprawdzić inżynierską metodą, jaka jest jakość obliczonego rozwiązania.

Najpierw musimy więc zdefiniować podstawowe obiekty wykorzystywane w procesie dyskretyzacji.

X = 2; % fizyczne rozmiary obszaru
k = 2.45; % stała dyfuzji

Nx = 199;
hx = X/(Nx+1);
px = hx*[0:Nx+1]'; 	% prawdziwe współrzędne węzłów dyskretyzacji
			% potrzebne m.in. do definicji prawej strony;
			% określone łącznie z brzegiem obszaru
Tx = lap1ddset(Nx,hx);	% macierz laplasjanu

% wyznaczamy funkcję prawej strony
m = 1; % punkt graniczny
pxl = px(px<m); pxr = px(px>=m); % węzły na lewo i prawo od "m"
Fx = [1./(2-pxl); 1./(pxr.^5)] / k;
plot(px,Fx); % sprawdzamy wizualnie poprawność Fx
Fx = Fx(2:end-1); % wybieramy tylko wartości w wewnętrznych węzłach

Zwróćmy uwagę na sposób generowania prawej strony: ponieważ musimy zdefiniować ją jedynie w wewnętrznych węzłach siatki, odwołujemy się jedynie do węzłów Fx(2:end-1). Dla większej skuteczności (i przejrzystości) kodu użyliśmy operatora potęgowania tablicowego.

Do zdefiniowania macierzy $T_{{N_{x}}}$ wykorzystaliśmy funkcję pomocniczą lap1ddset,

function T = lap1ddset(N, h)
	T = spdiags( repmat([-1 2 -1]./(h^2), N, 1), [-1,0,1], N, N);
end

Teraz wystarczy rozwiązać równanie i wizualizować wynik:

Ux = Tx \ Fx;
plot(px, [0; Ux; 0]); legend('U_{N_x}(x)'); grid on;

Pozostaje przekonać się, czy uzyskane przez nas rozwiązanie jest wystarczająco dokładne. W tym celu zastosujemy inżynierski sposób: sprawdzimy, jak bardzo różni się nasze rozwiązanie od rozwiązania uzyskanego na gęstszej siatce. Zapamiętamy więc Ux pod tymczasową nazwą,

TUx = Ux;

i ponownie uruchomimy symulację, startując jednak z $N_{x}=399$ (a więc, używając dwa razy drobniejszej siatki). Wtedy możemy ,,oszacować” popełniony błąd:

octave:12> norm(Ux(2:2:end-1)-TUx,Inf)
ans =  9.5256e-06

Zatem rozwiązanie uzyskane na dwukrotnie zagęszczonej siatce (a więc zapewne dokładniejsze, bo błąd aproksymacji maleje jak $O(h^{2})$ , por. wykład z numerycznych równań różniczkowych, rozdział dotyczący metody różnic skończonych dla równań eliptycznych) różni się w punktach starej siatki od oryginalnie uzyskanego jedynie o około $10^{{-5}}$ . To oczywiście tylko jest pewna sugestia co do błędu realnego, gdyż na pewno możemy powiedzieć jedynie, że w dyskretnej normie maksimum określonej w węzłach rzadszej siatki,

$\| u-U_{{N_{x}}}\| _{\infty}\leq\underbrace{\| U_{{N_{x}}}-U_{{2N_{x}+1}}\| _{\infty}}_{{\approx 10^{{-5}}\text{ (?)}}}+\underbrace{\| u-U_{{2N_{x}+1}}\| _{\infty}}_{{\text{wierzymy, że jest ,,małe'' wobec poprzedniego}}}.$

8.1.2.1. Testy

Nasza praca nie będzie zakończona, dopóki nie zweryfikujemy na kilku testowych przykładach, czy wykorzystana metoda produkuje rezultaty zgodne z oczekiwaniami. Minimalny zestaw testów mógłby tutaj wyglądać następująco:

Czy $T_{{N_{x}}}$ jest symetryczna? Czy jest trójdiagonalna tak, jak chcieliśmy?
Czy $h_{x}^{2}T_{{N_{x}}}$ jest macierzą o diagonali równej 2?
Ile wynosi $T_{{N_{x}}}U_{{N_{x}}}$ , gdy $U_{{N_{x}}}$ jest wektorem stałym? (Wynik powinien być zerowy z wyjątkiem pierwszego i ostatniego węzła.)
Czy jeśli $F_{{N_{x}}}=2$ , to rozwiązując zadanie dyskretne dostajemy w wyniku dyskretyzację funkcji $x(X-x)$ ? A dla $F_{{N_{x}}}=0$ ?
Funkcja $u(x)=\sin\left(\dfrac{2\pi}{X}x\right)$ spełnia równanie Poissona dla $f(x)=\left(\dfrac{2\pi}{X}\right)^{2}u(x)$ . Należy zbadać, czy błąd rozwiązań dyskretnych uzyskiwanych dla siatek o oczku $h_{x}$ , $h_{x}/2$ , $h_{x}/4$ , …, maleje zgodnie z teoretycznymi przewidywaniami dla tego zadania, por. wykład z numerycznych równań różniczkowych, rozdział o metodzie różnic skończonych dla równań eliptycznych.

Wykonanie tych testów pozostawiamy Czytelnikowi.

Ćwiczenie 8.3

Przeprowadź testy poprawności implementacji zarówno macierzy $T_{{N_{x}}}$ , jak i sposobu rozwiązywania jednowymiarowego równania Poissona.

Wskazówka: [pokaż]

Aby na przykład sprawdzić wynik działania operatora $T_{{N_{x}}}$ na wektorze stałym, wydamy komendę

Ux = ones(size(Fx));
s = Tx*Ux;
norm(s(2:end-1),Inf)
plot(px,[0;s;0]); grid on;

Obliczona powyżej norma fragmentu s powinna być równa zero, a wykres powinien mieć dwa piki w punktach tuż przy krańcach odcinka. Oczywiście, nie jest głupim pomysłem sprawdzenie wyniku dla kilku (np. małej i dużej) wartości $N_{x}$ .

8.1.3. Dyskretyzacja dwuwymiarowego laplasjanu

Chcąc rozwiązać dwuwymiarowe zadanie Poissona,

		$\displaystyle-\dfrac{\partial^{2}}{\partial x^{2}}u(x,y)-\dfrac{\partial^{2}}{\partial y^{2}}u(x,y)=f(x,y)\qquad\forall(x,y)\in(0,X)\times(0,Y)=\Omega,$		(8.10)
		$\displaystyle u(x,y)=0\qquad\forall(x,y)\in\partial\Omega,$		(8.11)

postąpimy analogicznie jak w poprzednim przypadku, korzystając z jednorodnej siatki węzłów na prostokącie $\Omega$ .

Jednorodna siatka na kwadracie, o identycznych kwadratowych oczkach.

Rys. 8.1. Regularna siatka na $\Omega$ .

Będziemy zatem mieli do czynienia z wartościami rozwiązania w punktach $(x_{i},y_{j})$ , gdzie

$x_{i}=i\cdot h_{x},\quad\text{dla }i=0,\ldots,N_{x}+1,\qquad\text{oraz}\qquad y_{j}=j\cdot h_{y},\quad\text{dla }j=0,\ldots,N_{y}+1,$

przy czym

$h_{x}=\dfrac{X}{N_{x}+1}\qquad\text{oraz}\qquad h_{y}=\dfrac{Y}{N_{y}+1}.$

Wartości rozwiązania dyskretnego $u_{{i,j}}$ będą odpowiadać wartościom rozwiązania w węłach dyskretyzacji $(x_{i},y_{j})$ , tzn. $u_{{i,j}}\approx u(x_{i},y_{j})$ .

Aby znaleźć przybliżone rozwiązanie (8.10), równanie różniczkowe zastąpimy odpowiadającym mu równaniem różnicowym, korzystając z tego samego co poprzednio sposobu aproksymacji drugiej pochodnej:

$\dfrac{\partial^{2}}{\partial x^{2}}u(x_{i},y_{j})\approx\dfrac{u_{{i-1,j}}-2u_{{i,j}}+u_{{i+1,j}}}{h_{x}^{2}}$

i analogicznie dla $\dfrac{\partial^{2}}{\partial y^{2}}(x_{i},y_{j})$ . Otrzymamy więc w końcu równanie różnicowe postaci

		$\displaystyle-\dfrac{u_{{i-1,j}}-2u_{{i,j}}+u_{{i+1,j}}}{h_{x}^{2}}-\dfrac{u_{{i,j-1}}-2u_{{i,j}}+u_{{i,j+1}}}{h_{y}^{2}}=f_{{i,j}}\qquad\forall i=1,\ldots,N_{x},\; j=1,\ldots,N_{y},$		(8.12)
		$\displaystyle u_{{0,j}}=u_{{N_{x}+1,j}}=u_{{i,0}}=u_{{i,N_{y}+1}}=0,\qquad\forall i=0,\ldots,N_{x}+1,\; j=0,\ldots,N_{y}+1,$		(8.13)

gdzie $f_{{i,j}}=f(x_{i},y_{j})$ . Mamy więc do wyznaczenia $N=N_{x}\cdot N_{y}$ niewiadomych, które w naturalny sposób układają się w macierz

$\begin{pmatrix}u_{{1,1}}&u_{{1,2}}&u_{{1,3}}&\ldots&u_{{1,N_{y}}}\\ u_{{2,1}}&u_{{2,2}}&u_{{2,3}}&\ldots&u_{{2,N_{y}}}\\ \vdots\\ u_{{N_{x},1}}&u_{{N_{x},2}}&u_{{N_{x},3}}&\ldots&u_{{N_{x},N_{y}}}\end{pmatrix}$

(zwróćmy uwagę na to, że ,,fizyczne” współrzędne są odwrócone: wartości odpowiadające kolejnym współrzędnym w kierunku osi $OX$ znajdują się w kolejnych wierszach jednej kolumny powyższej macierzy!). Podobnie jest z wartościami prawej strony równania. Aby jednak korzystać ze standardowych narzędzi Octave (na przykład operatora ,,\” rozwiązywania układów równań liniowych), nasze zadanie musimy sformułować w postaci, w której niewiadome (a także wektor prawej strony) są reprezentowane w postaci wektora, długości $N$ . Jakkolwiek można to robić na wiele sposobów, nam będzie najwygodniej uporządkować elementy kolumnami¹⁷Ze względu na to, że łatwo to osiągnąć w Octave; patrz niżej.:

$\displaystyle U_{N}=(\underbrace{u_{{1,1}},u_{{2,1}},u_{{3,1}},\ldots,u_{{N_{x},1}}}_{{\text{niewiadome pierwszej kolumny}}},\underbrace{u_{{1,2}},u_{{2,2}},u_{{3,2}},\ldots,u_{{N_{x},2}}}_{{\text{drugiej kolumny}}},\ldots,\underbrace{u_{{1,N_{y}}},u_{{2,N_{y}}},u_{{3,N_{y}}},\ldots,u_{{N_{x},N_{y}}}}_{{\text{ostatniej kolumny}}})^{T}.$

(8.14)

Analogicznie zdefiniowalibyśmy także wektor prawej strony $F_{N}$ . Będziemy więc szukać wektora $U_{N}$ , który spełnia liniowy układ równań o macierzowej postaci

$P_{N}\, U_{N}=F_{N},$

gdzie $P_{N}$ jest pięciodiagonalną macierzą o blokowej strukturze,

$P_{N}=\begin{pmatrix}T_{{N_{x}}}+\dfrac{2}{h_{y}^{2}}I&-\dfrac{1}{h_{y}^{2}}I&&&\\ -\dfrac{1}{h_{y}^{2}}I&T_{{N_{x}}}+\dfrac{2}{h_{y}^{2}}I&-\dfrac{1}{h_{y}^{2}}I&&\\ &\ddots&\ddots&\ddots&\\ &&-\dfrac{1}{h_{y}^{2}}I&T_{{N_{x}}}+\dfrac{2}{h_{y}^{2}}I&-\dfrac{1}{h_{y}^{2}}I\\ &&&-\dfrac{1}{h_{y}^{2}}I&T_{{N_{x}}}+\dfrac{2}{h_{y}^{2}}I\end{pmatrix}_{{N\times N}}.$

(8.15)

W kolejnych diagonalnych blokach $P_{N}$ znajdują się macierze trójdiagonalne $T_{{N_{x}}}+\dfrac{2}{h_{y}^{2}}I$ , gdzie $T_{{N_{x}}}$ jest znaną już nam macierzą dyskretyzacji jednowymiarowego laplasjanu. Jest tych bloków $N_{y}$ . Używając symbolu Kroneckera dla iloczynu tensorowego macierzy mamy po prostu

$P_{N}=I_{{N_{y}}}\otimes T_{{N_{x}}}+T_{{N_{y}}}\otimes I_{{N_{x}}}.$

Przypomnijmy na marginesie, że iloczyn tensorowy macierzy $A$ i $B$ jest macierzą o strukturze blokowej, postaci

$\begin{pmatrix}a_{{11}}B&a_{{12}}B&\ldots&a_{{1M}}B\\ a_{{21}}B&a_{{22}}B&\ldots&a_{{2M}}B\\ \vdots&\vdots&\ldots\\ a_{{N1}}B&a_{{N2}}B&\ldots&a_{{NM}}B\\ \end{pmatrix}.$

W Octave istnieje funkcja kron, która wyznacza iloczyn Kroneckera macierzy, przy czym zastosowana do macierzy rozrzedzonej daje w rezultacie także macierz rozrzedzoną. Dlatego definicja macierzy $P_{N}$ będzie w Octave wyjątkowo krótka:

hx = X/(Nx+1); hy = Y/(Ny+1);
Tx = lap1ddset(Nx, hx); Ix = speye(Nx, Nx);
Ty = lap1ddset(Ny, hy); Iy = speye(Ny, Ny);
P = kron(Iy, Tx) + kron(Ty, Ix);

Ćwiczenie 8.4

Zweryfikuj poprawność swojej implementacji macierzy $P_{N}$ .

8.1.4. Rozwiązanie równania Poissona w obszarze dwuwymiarowym

W przypadku macierzy $P_{N}$ , prosta strategia polegająca na eliminacji Gaussa (z wykorzystaniem dodatkowo faktu, że $P_{N}$ jest pasmowa, o szerokości pasma $2N_{x}$ ), dałaby nam szansę rozwiązać to równanie kosztem $O(N_{x}^{2}\, N_{y}^{2})$ , co jest, przy dużych wielkościach $N_{x}$ i $N_{y}$ , wartością trudną do zaakceptowania. Na szczęście, wbudowany Octave solver układów równań liniowych z macierzami rzadkimi, kryjący się za operatorem ,,\”, potrafi (czasem!) uporządkować niewiadome i równania tak, żeby wyraźnie obniżyć koszt rozwiązania.

Zatem, moglibyśmy w dalszym ciągu wyznaczać rozwiązanie poleceniem

U = P \ F;

o ile tylko $F_{N}$ w odpowiedni sposób przedstawimy w postaci długiego wektora F. Uzyskane rozwiązanie, U, będzie dla nas z pewnością bardziej użyteczne, jeśli ,,zwiniemy” je z powrotem do postaci macierzowej: U = reshape(U,Nx,Ny).

Jak to powinno wyglądać w szczegółach, jest treścią poniższego ćwiczenia.

Ćwiczenie 8.5

Wyznacz przybliżone rozwiązanie równania dyfuzji w środowisku anizotropowym,

		$\displaystyle\dfrac{\partial^{2}}{\partial x^{2}}u(x,y)-k\,\dfrac{\partial^{2}}{\partial y^{2}}u(x,y)=f(x,y),\qquad(x,y)\in(0,2)\times(0,3)=\Omega,$
		$\displaystyle u(x,y)=0,\qquad(x,y)\in\partial\Omega,$

gdzie $k=2.45$ oraz

$f(x,y)=\begin{cases}1,\qquad|x-1|\leq\frac{1}{2},|y-2|\leq\frac{1}{2},\\ 0,\qquad\text{w przeciwnym przypadku}.\end{cases}$

Narysuj wykres $u(x,y)$ i zbadaj następnie, jak rozwiązanie zmienia się wraz ze zmianą $k$ .

Rozwiązanie: [pokaż]

Zgodnie z wyżej opisanym schematem, musimy wyznaczyć macierz $P_{N}$ odpowiadającą równaniu różnicowemu

$-\dfrac{u_{{i-1,j}}-2u_{{i,j}}+u_{{i+1,j}}}{h_{x}^{2}}-k\,\dfrac{u_{{i,j-1}}-2u_{{i,j}}+u_{{i,j+1}}}{h_{y}^{2}}=f_{{i,j}}$

(zwróćmy uwagę na obecność współczynnika $k$ przy dyskretyzacji drugiej pochodnej względem $y$ ). Zatem, postępując zgodnie z opisaną metodą, po dyskretyzacji równania dostaniemy macierz

$P_{N}=I_{{N_{y}}}\otimes T_{{N_{x}}}+(k\, T_{{N_{y}}})\otimes I_{{N_{x}}}.$

Reszta jest już tylko (?) programowaniem. Musimy przy tym pamiętać, by we właściwym momencie z dwuwymiarowej tablicy wartości rozwiązania/prawej strony przechodzić do postaci rozwiniętej w wektor (lub na odwrót).

X = 2; Y = 3; k = 2.45;

Nx = 100; Ny = 200;
hx = X/(Nx+1); hy = Y/(Ny+1);

Tx = lap1ddset(Nx, hx); Ix = speye(Nx, Nx);
Ty = k*lap1ddset(Ny, hy); Iy = speye(Ny, Ny);
P = kron(Iy, Tx) + kron(Ty, Ix);

Chwila refleksji Jest jeszcze jeden ważny szczegół: funkcja prawej strony nie jest ciągła, zatem dokładne rozwiązanie nie może być klasy $C^{2}(\Omega)$ — nie będzie klasyczne. Owszem, możemy spodziewać się rozwiązań słabych — klasy $H^{1}_{0}(\Omega)$ — ale czy uprawnione w tym przypadku jest używanie aproksymacji różnicowej dla drugiej pochodnej, opartej na rozwinięciach Taylora?

Ten dylemat natury teoretycznej wydaje się być nieistotnym dla praktyka — w końcu przecież wystarczy użyć naszego sposobu i sprawdzić, ,,czy sobie poradzi”. Jednak kłopot polega na tym, że nie zawsze mamy wystarczające podstawy, by rozstrzygnąć, czy otrzymane rozwiązanie jest ,,sensowne”. Użycie kryterium czysto estetycznego (,,ładnie wygląda”) niestety często prowadzi do pochopnej akceptacji całkowicie fałszywych wyników.

Zatem wstępne zastanowienie się nad jakością używanych narzędzi aproksymacji problemu nie jest bez znaczenia, bo użycie błędnej/niewłaściwej dyskretyzacji może w efekcie doprowadzić nas do efektownych wizualizacji, które nie będą miały wiele wspólnego z realnym rozwiązaniem. Zachęcając Czytelnika do poszukania odpowiedzi¹⁸Okazuje się, że nasza dyskretyzacja ma sens nawet w tym przypadku: gdy potraktować ją jako dyskretyzację metodą elementu skończonego z nienajlepszą (ale wciąż sensowną) aproksymacją wektora prawej strony. Następnym etapem rozwiązywania tego zadania byłoby więc ,,podrasowanie” wektora prawej strony. na to pytanie na przykład w naszym ulubionym rozdziale o metodzie różnic skończonych dla równań eliptycznych w wykładzie z numerycznych równań różniczkowych, teraz — na własne ryzyko — (a przy tym dla lepszego rozeznania) sprawdzimy jednak po prostu, ,,co wyjdzie”.

Na początek wygodnie będzie nam utworzyć macierz $F=\begin{pmatrix}F_{{i.j}}\end{pmatrix}$ określoną we wszystkich węzłach dyskretyzacji (łącznie z brzegowymi):

px = hx*[0:Nx+1]'; 	% prawdziwe współrzędne węzłów dyskretyzacji
py = hy*[0:Ny+1]';
[pX, pY] = meshgrid(px,py); % tablica węzłów do manipulacji tablicowej

F = zeros(size(pX));
[i,j] = find( (abs(pX-1) <= 0.5) & (abs(pY-2) <= 0.5) );
F(i,j) = 1;
F = permute(F, [2 1]); % zamienia ze sobą wiersze i kolumny

Ostatni krok może wydać się nam niejasny: dlaczego musimy zamienić ze sobą wiersze i kolumny macierzy $F$ ? Przyczyną jest struktura, jaką mają macierze pX, pY, które Octave generuje poleceniem meshgrid. Aby ułatwić wizualizację takich macierzy, indeksy odpowiadające zmiennej $x$ odpowiadają tam numerom kolumn — a nie, tak, jak my sobie to założyliśmy, wierszy. Określając więc wartości F na podstawie pX i pY, dostajemy F o zamienionych ze sobą wierszach i kolumnach. Dlatego na koniec musimy to odwrócić poleceniem permute¹⁹Tak, moglibyśmy też użyć transpozycji, ale nie uogólniłoby się to na przypadek trójwymiarowy, którym za niedługo będziemy się zajmować..

Chcąc teraz wizualizować macierz F, musimy ponownie zamienić ze sobą wymiary: polecenie imagesc oczekuje macierzy wartości w węzłach takiej, jak np. produkowana poleceniem meshgrid.

imagesc(px, py, permute(F, [2 1]));
axis('xy','square');
title('Prawa strona');

Na zakończenie musimy z tablicy F usunąć niepotrzebne wartości odpowiadające węzłom leżącym na brzegu, a następnie ,,rozprostować” F do wektora:

F = F(2:Nx+1,2:Ny+1);
F = F(:);

Ostatnie polecenie tworzy wektor, którego elementy są kolejnymi elementami kolejnych kolumn macierzy F — dokładnie tak, jak umawialiśmy się w (8.14).

Teraz wystarczy zadanie rozwiązać, a wynik na powrót zwinąć do tablicy i wizualizować.

U = zeros(Nx+2,Ny+2);
U(2:Nx+1, 2:Ny+1) = reshape(P \ F, Nx, Ny);
imagesc(px, py, permute(U, [2 1]));
axis('xy','square');
title('Rozwiazanie');

Na zakończenie powinniśmy ocenić (w analogiczny sposób jak w przypadku jednowymiarowym), jakiego rzędu błąd popełniliśmy stosując naszą dyskretyzację. Wykonanie odpowiednich testów pozostawiamy Czytelnikowi.

8.1.5. Dyskretyzacja i rozwiązanie równania Poissona w obszarze trójwymiarowym

Czytelnik zechce sam sprawdzić, że dyskretyzacja operatora Laplace'a w kostce $\Omega=(0,X)\times(0,Y)\times(0,Z)$ różnicami skończonymi na regularnej siatce o oczku $(h_{x},h_{y},h_{z})$ i $N_{x}\times N_{y}\times N_{z}$ wewnętrznych węzłach, prowadzi do zadania różnicowego

$\displaystyle-\dfrac{u_{{i-1,j,k}}-2u_{{i,j,k}}+u_{{i+1,j,k}}}{h_{x}^{2}}-\dfrac{u_{{i,j-1,k}}-2u_{{i,j,k}}+u_{{i,j+1,k}}}{h_{y}^{2}}-\dfrac{u_{{i,j,k-1}}-2u_{{i,j,k}}+u_{{i,j,k+1}}}{h_{z}^{2}}=f_{{i,j,k}}$

(8.16)

dla $i=1,\ldots,N_{x},\; j=1,\ldots,N_{y},\; k=1,\ldots,N_{z}$ , z odpowiednimi warunkami brzegowymi. Oczywiście $u_{{i,j,k}}\approx u(ih_{x},jh_{y},kh_{z})$ . Mamy więc do wyznaczenia $N=N_{x}\cdot N_{y}\cdot N_{z}$ niewiadomych wewnątrz obszaru, które w naturalny sposób układają się w trójwymiarową tablicę

$\begin{pmatrix}u_{{i,j,k}}\end{pmatrix}_{{N_{x}\times N_{y}\times N_{z}}}.$

Aby wciąż korzystać ze standardowych narzędzi Octave jak w przypadkach dwu- i jednowymiarowym, nasze zadanie musimy sformułować w postaci, w której niewiadome tworzą jeden wektor długości $N$ . Z objaśnionych wcześniej powodów, zrobimy to ,,kolumnami” — to znaczy tak, by w uzyskanym wektorze najszybciej zmieniała się pierwsza współrzędna, potem druga, a najwolniej — trzecia:

	$\displaystyle U_{N}=($	$\displaystyle\underbrace{u_{{1,1}},u_{{2,1,1}},u_{{3,1,1}},\ldots,u_{{N_{x},1,1}}}_{{\text{niewiadome pierwszej kolumny i pierwszego plastra}}},\underbrace{u_{{1,2,1}},u_{{2,2,1}},u_{{3,2,1}},\ldots,u_{{N_{x},2,1}}}_{{\text{drugiej kolumny i pierwszego plastra}}},\ldots,$
		$\displaystyle\underbrace{u_{{1,N_{y},1}},u_{{2,N_{y},1}},u_{{3,N_{y},1}},\ldots,u_{{N_{x},N_{y},1}}}_{{\text{ostatniej kolumny i pierwszego plastra}}},\ldots,$
		$\displaystyle\ldots,$
		$\displaystyle\underbrace{u_{{1,N_{y},N_{z}}},u_{{2,N_{y},N_{z}}},u_{{3,N_{y},N_{z}}},\ldots,u_{{N_{x},N_{y},N_{z}}}}_{{\text{ostatniej kolumny i ostatniego plastra}}})^{T}.$

Będziemy szukać wektora $U_{N}$ , który spełnia liniowy układ równań

$S_{N}\, U_{N}=F_{N},$

gdzie $S_{N}$ jest siedmiodiagonalną macierzą o blokowej strukturze,

$S_{N}=I_{{N_{z}}}\otimes I_{{N_{y}}}\otimes T_{{N_{x}}}+I_{{N_{z}}}\otimes T_{{N_{y}}}\otimes I_{{N_{x}}}+T_{{N_{z}}}\otimes I_{{N_{y}}}\otimes I_{{N_{x}}}.$

(8.17)

Implementacja takiej macierzy, a nastepnie jej rozwiązanie za pomocą standardowego bezpośredniego solvera z Octave przebiega w sposób w pełni analogiczny do przypadku dwuwymiarowego. Niestety, gdy ledwie dziesięciokrotnie zmniejszymy rozmiar oczka siatki w każdym kierunku, liczba niewiadomych rośnie aż tysiąckrotnie. Otrzymane zadanie będzie więc bardzo duże: faktycznie, już dla umiarkowanego $N_{x}=N_{y}=N_{z}=100$ dochodzimy do $N=10^{6}$ (!) Dodatkowo, struktura naszej (bardzo rozrzedzonej) macierzy jest na tyle skomplikowana dla bezpośredniego solvera opartego na eliminacji Gaussa, że szybko staje się on nieskuteczny.

Ćwiczenie 8.6

Rozwiąż zagadnienie Poissona

	$\displaystyle-\Delta u$	$\displaystyle=f\qquad\text{w }\Omega=(0,1)\times(0,2)\times(0,3),$
	$\displaystyle u$	$\displaystyle=0\qquad\text{na }\partial\Omega,$

gdzie $f$ przyjmuje wartość 1, zaburzoną o losowe wartości z przedziału $[-0.1,0.1]$ . Wizualizuj rozwiązanie.

Rozwiązanie: [pokaż]

Tak sformułowane zadanie ma charakter czysto jakościowy, a nie ilościowy — rzeczywiście, skoro $f$ jest poddana losowemu zaburzeniu, to trudno spodziewać się, by błąd popełniany z powodu dyskretyzacji miał inny niż losowy charakter i dlatego jakikolwiek konkretny wynik sam w sobie jest mało użyteczny: wszak rodzielczość użytej przez nas siatki nigdy nie będzie wystarczająca, by uchwycić wszelkie wahnięcia $f$ ! Jednak zadanie będzie dla nas dobrym pretekstem do tego, by

pokazać wygładzający efekt dyfuzji
zobaczyć ograniczenia bezpośredniego solvera układów równań z macierzami rzadkimi na jakimkolwiek konkretnym przykładzie.

X = 1; Y = 2; Z = 3; % fizyczne rozmiary kostki
Nx = 30; Ny = 30; Nz = 30;

hx = X/(Nx+1);
hy = Y/(Ny+1);
hz = Z/(Nz+1);

px = hx*[0:Nx+1]; py = hy*[0:Ny+1]; pz = hz*[0:Nz+1]; % prawdziwe współrzędne węzłów (lacznie z brzegiem obszaru)

% macierze pomocnicze dla zadan jednowymiarowych
Tx = lap1ddset(Nx,hx); Ty = lap1ddset(Ny,hy); Tz = lap1ddset(Nz,hz);
Ix = speye(size(Tx)); Iy = speye(size(Ty)); Iz = speye(size(Tz));

% macierz laplasjanu
S = kron(Iz, kron(Iy,Tx)) + kron(Iz, kron(Ty, Ix)) + kron(Tz, kron(Iy, Ix));

% prawa strona i jej wizualizacja
F = 1 + 0.1*(2*rand(Nx+2,Ny+2,Nz+2)-1);
slice(px, py, pz, permute(F, [2 1 3]), X/2,Y/2, Z/2);
xlabel('x'); ylabel('y'); zlabel('z'); grid on; title('Prawa strona');
pause(1);

F = F(2:Nx+1,2:Ny+1,2:Nz+1); % ograniczamy się do węzłów wewnętrznych

% wyznaczamy rozwiązanie
U = zeros(Nx+2,Ny+2,Nz+2);
U(2:Nx+1,2:Ny+1,2:Nz+1) = reshape(S \ F(:),Nx,Ny,Nz);

% wizualizacja rozwiązania
slice(px, py, pz, permute(U, [2 1 3]), X/2,Y/2, Z/2);
xlabel('x'); ylabel('y'); zlabel('z'); grid on; title('Rozwiazanie');
pause(1);

Zwróćmy uwagę na to, jak grubą siatkę wzięliśmy w naszym przykładowym kodzie. Gdyby bowiem zagęścić ją dwukrotnie, czas potrzebny na uzyskanie rozwiązania znacząco wydłuża się; próba wyznaczenia rozwiązania na siatce $80\times 80\times 80$ kończy się²⁰Na komputerze z 2GB pamięci RAM. Na szybkiej maszynie z ,,nieograniczonymi” zasobami (wieloprocesorowy Intel Xeon 2.4GHz, 64GB RAM), rozwiązanie udaje się policzyć, owszem, wykorzystując w tym celu (w MATLABie) ponad 8.5 GB pamięci i 400 sekund czasu rzeczywistego. komunikatem o błędzie z powodu braku wystarczających zasobów.

8.1.5.1. Krytyka standardowego solvera bezpośredniego

Jakkolwiek stosowany przez nas solver bezpośredni całkiem dobrze sobie radzi z zadaniami dwuwymiarowymi, w których liczba węzłów sięga nawet $500\times 500$ , to już dla zadań rozmiaru $1000\times 1000$ staje się mało efektywny, a jego wymagania pamięciowe stają się duże. Sytuacja staje się jeszcze bardziej dramatyczna, gdy przejdziemy do zadań trójwymiarowych, co mieliśmy okazję oglądać w ostatnim przykładzie. Jasne jest, że stosowane przez bezpośredni solver algorytmy reorderingu niewiadomych i równań są nieskuteczne w przypadku dyskretyzacji trójwymiarowych — nawet tak regularnych jak nasza.

8.1.5.2. Solver oparty na FFT

W przypadku tak specyficznej macierzy, jaką jest macierz laplasjanu, można na szczęście wskazać specjalne metody jej rozwiązywania. To typowa i charakterystyczna dla obliczeń naukowych sytuacja, kiedy do konkretnego zadania można dobrać niestandardową metodę, która będzie znacznie skuteczniejsza od metod ogólnego przeznaczenia. Często taka metoda będzie motywowana fizycznymi chechami zadania; w naszym przypadku lepsza metoda będzie wynikać wyłącznie z matematycznych własności zadania.

Otóż okazuje się, że znane są jawne wzory na czynniki rozkładu macierzy $T_{N}$ na wartości i wektory własne,

$T_{N}=Q_{N}\Lambda _{N}Q_{N}^{T},$

skąd wynika, że

w przypadku dwuwymiarowym:

$P_{N}=Q_{{N_{y}}}Q_{{N_{x}}}(\underbrace{I_{{N_{y}}}\otimes\Lambda _{{N_{x}}}+\Lambda _{{N_{y}}}\otimes I_{{N_{x}}}}_{{\text{macierz diagonalna}}})Q_{{N_{x}}}^{T}Q_{{N_{y}}}^{T}.$
w przypadku trójwymiarowym:

$S_{N}=Q_{{N_{z}}}Q_{{N_{y}}}Q_{{N_{x}}}(\underbrace{I_{{N_{z}}}\otimes I_{{N_{y}}}\otimes\Lambda _{{N_{x}}}+I_{{N_{z}}}\otimes\Lambda _{{N_{y}}}\otimes I_{{N_{x}}}+\Lambda _{{N_{z}}}\otimes I_{{N_{y}}}\otimes I_{{N_{x}}}}_{{\text{macierz diagonalna}}})Q_{{N_{x}}}^{T}Q_{{N_{y}}}^{T}Q_{{N_{z}}}^{T}.$

Rozwiązanie układu z macierzą diagonalną kosztuje tylko $N$ operacji, natomiast mnożenia przez macierze $Q_{N}$ można łatwo zrealizować (niskim kosztem) za pomocą algorytmu FFT, dostępnego w Octave. Więcej informacji na ten temat (w tym — informację o jeszcze szybszym algorytmie!) można znaleźć w wykładach z Matematyki Obliczeniowej II.

Przykład 8.1 (Trójwymiarowe równanie Poissona przez FFT)

Pokażemy, jak rozwiązywać trójwymiarowe równanie Poissona (por. ćwiczenie 8.6) z wykorzystaniem solvera opartego na FFT. Jest to niewielka modyfikacja kodu z ćwiczenia 8.6: niepotrzebne fragmenty zostały wykomentowane, a nowo wstawiane części oznakowano komentarzem NOWE.

X = 1; Y = 2; Z = 3; % fizyczne rozmiary kostki
Nx = 100; Ny = 100; Nz = 100; scale = 0.1;

hx = X/(Nx+1);
hy = Y/(Ny+1);
hz = Z/(Nz+1);

px = hx*[0:Nx+1]; py = hy*[0:Ny+1]; pz = hz*[0:Nz+1]; % prawdziwe wspolrzedne wezlow (lacznie z brzegiem obszaru)

% macierze pomocnicze dla zadan jednowymiarowych
% Tx = lap1ddset(Nx,hx); Ty = lap1ddset(Ny,hy); Tz = lap1ddset(Nz,hz); % potrzebne tylko dla metody bezposredniej
Ix = speye(Nx,Nx); Iy = speye(Ny,Ny); Iz = speye(Nz,Nz);

% NOWE!
Lx = spdiags(fftsetd(Nx,hx),0,Nx,Nx); Ly = spdiags(fftsetd(Ny,hy),0,Ny,Ny); Lz = spdiags(fftsetd(Nz,hz),0,Nz,Nz); % macierze diagonalne, potrzebne tylko dla metody FFT

% macierz laplasjanu, tylko dla metody bezposredniej
% S = kron(Iz, kron(Iy,Tx)) + kron(Iz, kron(Ty, Ix)) + kron(Tz, kron(Iy, Ix));

% NOWE!
% macierz diagonalna zadania 3D, tylko dla metody FFT
global Lambda; % gloablna, bo wykorzystywana przez fftsolve3dd()
Lambda = kron(Iz, kron(Iy,Lx)) + kron(Iz, kron(Ly, Ix)) + kron(Lz, kron(Iy, Ix));

% prawa strona
F = 1 + scale*(2*rand(Nx+2,Ny+2,Nz+2)-1);
slice(px, py, pz, permute(F, [2 1 3]), X/2,Y/2, Z/2);
xlabel('x'); ylabel('y'); zlabel('z'); grid on; title('Prawa strona');
pause(1);

F = F(2:Nx+1,2:Ny+1,2:Nz+1);

% tylko dla metody bezposredniej
% U = zeros(Nx+2,Ny+2,Nz+2);
% U(2:Nx+1,2:Ny+1,2:Nz+1) = reshape(S \ F(:),Nx,Ny,Nz); % rozwiazanie

% NOWE!
% macierz diagonalna zadania 3D, tylko dla metody FFT
U = zeros(Nx+2,Ny+2,Nz+2);
tic; U(2:Nx+1,2:Ny+1,2:Nz+1) = fftsolve3dd(F); toc % rozwiazanie FFT

% wizualizacja
slice(px, py, pz, permute(U, [2 1 3]), X/2,Y/2, Z/2);
xlabel('x'); ylabel('y'); zlabel('z'); grid on; title('Rozwiazanie');
pause(1);

Powyższy skrypt korzysta z pewnych dodatkowych funkcji: fftsolve3dd i fftsetd; poniżej ich listingi.

function V = fftsolve3dd(F)
% Rozwiazuje, korzystajac z FFT, uklad rownan liniowych Lap*V = F, 
% gdzie V,F sa macierzami wartosci na siatce kwadratowej.
% Warunki brzegowe Dirichleta: jednorodne.
global Lambda; % wartosci wlasne

[Nx, Ny, Nz] = size(F);

F = fftmult3d(F);	%F = ZX*F*ZY;

V = Lambda\F(:);

V = reshape(V, Nx, Ny, Nz);

V = fftmult3d(V);	%V = ZX*V*ZY;

end

function Lambda = fftsetd(N,h)
%
% Wartosci wlasne dla 1-d Laplasjanu z war brzeg. Dirichleta

Lambda = (4/h^2)*(sin(0.5*pi*[1:N]'/(N+1))).^2;	% exact eigenvalues
end

Nagle… realne staje się szybkie wyznaczenie rozwiązania zadania trójwymiarowego, nie tylko na siatce $100\times 100\times 100$ (to milion niewiadomych!) w czasie poniżej 1 sekundy (sic!) na komputerze Core 2 Duo 2.4 GHz z pamięcią 2GB. Nawet zadanie postawione na siatce $200\times 200\times 200$ wciąż liczy się w więcej niż przyzwoitym czasie: poniżej 7 sekund!

Kluczowa dla szybkości działania solvera jest funkcja fftsolve3dd, wykorzystująca umiejętność Octave przykładania FFT ,,wzdłuż” zadanego wymiaru macierzy.

Ćwiczenie 8.7

Zmodyfikuj powyższy kod oraz funkcję fftsolve3dd tak, by mogła działać w przypadku dyskretyzacji dwuwymiarowego równania Poissona.

Quiz 8.1

Czy warto zastosować powyższy algorytm do rozwiązywania jednowymiarowego zadania Poissona?
Tak

Nie
Dobrze Metoda oparta na eliminacji Gaussa, wykorzystująca strukturę macierzy trójdiagonalnej, jest w tym wypadku optymalna
Źle Pomyśl, ile kosztuje rozwiązanie układu równań z macierzą trójdiagonalną
Czy w przypadku dwuwymiarowym opisana powyżej metoda oparta na FFT będzie dawała rozwiązanie obarczone większym błędem niż w sytuacji gdy zastosujemy w Octave standardowy operator rozwiązywania równań liniowych?
Tak

Nie
Dobrze Odpowiedź, przynamniej teoretycznie, brzmi: nie. Obie metody rozwiązują ten sam układ równań i obie są metodami bezpośrednimi: w arytmetyce dokładnej dają dokładne rozwiązanie. Ponieważ w praktyce obliczenia wykonujemy w arytmetyce podwójnej precyzji, a każda z metod wykonuje inne rachunki, należy spodziewać się niewielkiej różnicy, ,,na poziomie precyzji arytmetyki”. Większy wpływ na wynik może mieć złe uwarunkowanie samej macierzy, a nie różnice między solverami. Patrz także [11]
Źle To jeszcze nie uzasadnienie, ale… która z metod wyznaczy rozwiązanie niższym kosztem. Zajrzyj do [11]

[sprawdź odpowiedzi]

Ćwiczenie 8.8

Sprawdź, czy metoda iteracyjna — taka jak na przykład PCG, o której więcej dowiesz się z wykładu z Matematyki Obliczeniowej II, a która jest implementowana w Octave komendą pcg — jest rozsądną alternatywą dla solvera FFT.

Wskazówka: [pokaż]

W przypadku zadania trójwymiarowego należy użyć komendy

U = zeros(Nx+2,Ny+2,Nz+2);
[Uin, flag, relres, iter, resvec, eigest] = ...
	pcg(S,F(:), (1e-3)*min([hx,hy,hz])^2, 50);
U(2:Nx+1,2:Ny+1,2:Nz+1) = reshape(Uin,Nx,Ny,Nz);
disp(['Uwarunkowanie S = ', num2str(eigest(2)/eigest(1))]);

Jakie jest uwarunkowanie tej macierzy? (Odpowiedź na to pytanie znajdziesz na przykład znów w wykładzie z Matematyki Obliczeniowej II.) Co to oznacza dla solvera iteracyjnego, takiego jak PCG?

8.2. Równanie dyfuzji

Wykorzystywany przez nas solver FFT ma niestety poważne ograniczenia, które uniemożliwiają zastosowanie go wprost do zadania (8.1), w którym mamy co prawda do czynienia z jednorodną siatką na kostce $\Omega$ , ale jednocześnie współczynnik dyfuzji nie jest stały w całym obszarze: aby go zastosować, musieliśmy założyć, że mamy do czynienia z równaniem Poissona, bo tylko dla niego znaliśmy wartości i wektory własne odpowiednich macierzy.

Okazuje się wszakże, iż możemy efektywnie skorzystać z FFT nawet w przypadku zmieniającego się współczynnika dyfuzji. Jednak wtedy musimy podejść do zadania z jeszcze innej strony: używając metody iteracyjnej PCG do rozwiązywania zadania i wykorzystując jednocześnie bezpośredni solver oparty na FFT w charakterze wydajnego operatora ściskającego²¹Więcej szczegółów na temat PCG, ściskania i macierzy ściskających można poznać na wykładzie z Matematyki Obliczeniowej II., dzięki czemu metoda PCG będzie szybko zbieżna!

Zwróćmy uwagę, że wiele metod iteracyjnych — w tym metoda PCG — rozwiązywania układu równań liniowych

$Ax=b$

nie wymaga wcale konstrukcji macierzy $A$ ; w zupełności wystarcza im funkcja, wyznaczająca dla zadanego wektora $v$ ilocznyn $Av$ (takie metody iteracyjne nazywamy metodami operatorowymi). W naszym przypadku będzie to wielkim ułatwieniem, gdyż obliczenie działania dyskretnego operatora dyfuzji na zadanej funkcji siatkowej będzie łatwe do osiągnięcia — w przeciwieństwie do jawnej konstrukcji macierzy takiego operatora. To jeszcze jeden powód, by zastosować metodę PCG.

Aby przetestować nasze nowe podejście, zaczniemy od przypadku jednowymiarowego, uwzględniając (wreszcie!) przypadek, gdy warunek brzegowy nie jest jednorodny.

8.2.1. Dyskretyzacja zadania jednowymiarowego

Gdy $\Omega$ jest odcinkiem $(0,X)$ , (8.1) staje się równaniem różniczkowym zwyczajnym z warunkiem brzegowym:

		$\displaystyle-\dfrac{d}{dx}(a(x)\,\dfrac{d}{dx}u(x))=f(x)\qquad\forall x\in(0,X),$		(8.18)
		$\displaystyle u(0)=g(0),\quad u(X)=g(X).$		(8.19)

Aby znaleźć jego przybliżone rozwiązanie, możemy na przykład zastąpić równanie różniczkowe odpowiadającym mu równaniem różnicowym. Okazuje się, że prostoduszne

$-\dfrac{a_{{i-1}}u_{{i-1}}-2a_{{i}}u_{i}+a_{{i+1}}u_{{i+1}}}{h_{x}^{2}}=f_{i}$

dość kiepsko aproksymuje nasze równanie, a znacznie lepsze własności (wiedza nie szkodzi!) ma aproksymacja postaci

		$\displaystyle-\dfrac{a_{{i-1/2}}u_{{i-1}}-(a_{{i-1/2}}+a_{{i+1/2}})u_{i}+a_{{i+1/2}}u_{{i+1}}}{h_{x}^{2}}=f_{i}\qquad\forall i=1,\ldots,N_{x},$		(8.20)
		$\displaystyle u_{0}=g_{0},\qquad u_{{N_{x}+1}}=g_{{N_{x}+1}},$		(8.21)

gdzie za $a_{{i+1/2}}$ można wziąć na przykład

$a_{{i+1/2}}=\dfrac{a_{{i}}+a_{{i+1}}}{2}\qquad\text{lub}\qquad a_{{i+1/2}}=a(x_{{i+1/2}}).$

gdzie, zgodnie z intuicją, $x_{{i+1/2}}=\dfrac{x_{{i}}+x_{{i+1}}}{2}$ . Widzimy więc, że obliczenie działania dyskretnego operatora dyfuzji $D_{{N_{x}}}:\mathbb{R}^{{N_{x}+2}}\rightarrow\mathbb{R}^{{N_{x}}}$ na wektorze $U_{{N_{x}}}$ (dane wzorem (8.20)), wymaga znajomości wartości także $U_{{N_{x}}}$ na brzegu obszaru i prowadzi do wyznaczenia wartości wyniku jedynie w węzłach wewnętrznych.

8.2.1.1. Uwzględnienie warunków brzegowych

W naszym zadaniu poszukujemy rozwiązania postaci

$U_{{N_{x}}}=(g_{0},\underbrace{u_{1},\ldots,u_{{N_{x}}}}_{{\text{realne niewiadome}}},g_{{N_{x}+1}})^{T}$

takiego, że

$D_{{N_{x}}}U_{{N_{x}}}=(f_{1},\ldots,f_{{N_{x}}})^{T}.$

Powyższe zadanie, z niejednorodnym warunkiem brzegowym, bardzo łatwo sprowadzić do zadania z warunkiem jednorodnym. Rzeczywiście, oznaczając

$\tilde{U}_{{N_{x}}}=U_{{N_{x}}}-(g_{0},\underbrace{0,\ldots,0}_{{\text{wewnętrzne}}},g_{{N+1}})^{T}$

mamy, że $\tilde{U}_{{N_{x}}}$ na brzegu przyjmuje wartości równe zero oraz spełnia

$\displaystyle D_{{N_{x}}}\tilde{U}_{{N_{x}}}=(f_{1},\ldots,f_{{N_{x}}})^{T}-D_{{N_{x}}}(g_{0},0,\ldots,0,g_{{N+1}})^{T}.$

(8.22)

W ten sposób wystarczy zmodyfikować prawą stronę równania, by potem przez rozwiązanie równania z jednorodnym warunkiem brzegowym wyznaczyć wartości rozwiązania we wszystkich punktach wewnętrznych siatki.

Oto więc, jak można byłoby zaimplementować działanie operatora $D_{{N_{x}}}$ w Octave:

function F = multdiff1dd(Uin, h, A, U)
%
% wyznacza F = -(a(x) u'(x))', dla dyskretyzacji
%
% Uin - wartości u(x) we wnętrzu siatki
% h - długość kroku siatki
% A - wartości a(x) we wszystkich węzłach pośrednich siatki
%	domyślnie a(x) = 1
% U - wartości zerowego rozszerzenia warunku brzegowego na wszystkie węzły siatki
%	domyślnie U(x) = 0

if nargin < 4 % domyślnie jednorodny warunek brzegowy
	U = zeros(length(Uin)+2,1);
end

[Nx, Ny] = size(U);

if nargin < 3 % domyślnie a(x) = 1
	Ax = ones(Nx-1,1);
else
	Ax = A;
end

ix = 2:Nx-1; % indeksy odpowiadajace wartosciom z wnetrza

U(2:Nx-1) = Uin;
F = (-U(ix-1).*Ax(ix-1) ...
	+ U(ix).*(Ax(ix-1)+Ax(ix)) ...
	- U(ix+1).*Ax(ix))/(h(1)^2);
end

Jeśli więc chcemy wyznaczyć wynik działania $D_{{N_{x}}}$ na

tablicy U o wartościach zadanych wyłącznie w wewnętrznych węzłach siatki i domyślnej wartości zero na brzegu, użylibyśmy wywołania
```
F = multdiff1dd(U, hx, Ax);
```
tablicy o wartościach w wewnętrznych węzłach siatki równej zero i wartościach na brzegu określonych przez tablicę G (rozmiaru $N_{x}+2$ ), użylibyśmy wywołania
```
F = multdiff1dd(zeros(Nx,1), hx, Ax, G);
```

Ostatnie wywołanie pasuje jak ulał do wyznaczenia poprawki na prawą stronę $D_{{N_{x}}}(g_{0},0,\ldots,0,g_{{N+1}})^{T}$ , por. (8.22).

Ćwiczenie 8.9

Pamiętając o konieczności nieustannego weryfikowania poprawności zaimplementowanych funkcji, przygotuj i przeprowadź serię testów sprawdzających funkcję multdiff1dd.

Wskazówka: [pokaż]

8.2.2. Rozwiązanie równania dyfuzji w obszarze jednowymiarowym

Zgodnie z wcześniej przyjętą strategią, zastosujemy iteracyjną metodę rozwiązywania równania. Dla ustalenia uwagi przyjmieny, że $a(x)=x+1$ oraz $f(x)=g(x)=1$ .

A = pX+1; % definicja wartości współczynnika dyfuzji w węzłach siatki
Ax = 0.5*(A(1:end-1) + A(2:end));
multdiff = @(x) multdiff1dd(x, hx, Ax); % funkcja przykładająca operator dyfuzji z zerowym warunkiem brzegowym
precond = @(x) Tx\x; % jako operator ściskający wybieramy rozwiązanie r-nia Poissona

F = ones(Nx,1);
U = ones(Nx+2,1);
F = F - multdiff1dd(zeros(Nx,1), hx, Ax, U); % uwzględnienie warunku brzegowego
[Uin, flag, relres, iter, resvec, eigest] = ...
		pcg(multdiff, F(:), (1e-3)*hx^2, 50, precond);
U(2:Nx+1) = Uin;

Ćwiczenie 8.10

Czy w przypadku jednowymiarowej dyfuzji nie można byłoby w prosty sposób użyć solvera bezpośredniego?

Rozwiązanie: [pokaż]

Ćwiczenie 8.11

Pamiętając o konieczności nieustannego weryfikowania poprawności zaimplementowanych funkcji, przygotuj i przeprowadź serię testów sprawdzających opracowaną metodę.

Wskazówka: [pokaż]

Ćwiczenie 8.12

Funkcja mnożąca wektor przez operator ściskający powinna być możliwie tania, gdyż będziemy ją wywoływać wiele razy. Nasza funkcja:

precond = @(x) Tx \ x;

ma tę wadę, że za każdym wywołaniem dokonuje rozkładu LU macierzy $T_{x}$ . Pomyśl, jak to usprawnić.

Rozwiązanie: [pokaż]

Wyznaczając na etapie preprocessingu rozkład LU (lub Cholesky'ego)

Rx = spchol(Tx); % rozkład Cholesky'ego: Rx'*Rx = Tx
precond = @(x) Rx \ (Rx' \ x);

8.2.3. Dyskretyzacja dwuwymiarowego operatora dyfuzji

Gdybyśmy od razu chcieli atakować zadanie (8.1) w przypadku dwu- lub więcejwymiarowym, mogłoby nam być dosyć trudno. Teraz jednak, po przepracowaniu przypadku jednowymiarowego, sprawa powinna potoczyć się wartko. Będziemy musieli pamiętać o dwoistości spojrzenia na niewiadome: raz będziemy je traktować jako (dwuwymiarową!) tablicę wartości w węzłach (tak będzie nam wygodniej zapisać działanie operatora dyfuzji), innym razem będziemy musieli rozwinąć tablicę w jeden długi wektor (tak działa funkcja pcg w Octave).

Przez analogię do przypadku jednowymiarowego, nasza dyskretyzacja będzie miała postać:

		$\displaystyle-\dfrac{a_{{i-1/2,j}}u_{{i-1,j}}-(a_{{i-1/2,j}}+a_{{i+1/2,j}})u_{{i,j}}+a_{{i+1/2,j}}u_{{i+1,j}}}{h_{x}^{2}}$
		$\displaystyle\qquad-\dfrac{a_{{i,j-1/2}}u_{{i,j-1}}-(a_{{i,j-1/2}}+a_{{i,j+1/2}})u_{{i,j}}+a_{{i,j+1/2}}u_{{i,j+1}}}{h_{y}^{2}}$
		$\displaystyle\qquad=f_{{i,j}}\qquad\forall i=1,\ldots,N_{x},\; j=1,\ldots,N_{y},$

z warunkiem brzegowym

$u_{{0,j}}=g_{{0,j}},\quad u_{{N_{x}+1,j}}=g_{{N_{x}+1,j}},\quad u_{{i,0}}=g_{{i,0}},\quad u_{{i,N_{y}+1}}=g_{{i,N_{y}+1}},\qquad\forall i=1,\ldots,N_{x},\; j=1,\ldots,N_{y}.$

Mamy do wyznaczenia $N=N_{x}\cdot N_{y}$ niewiadomych, które w naturalny sposób układają się w macierz

$U_{N}=\begin{pmatrix}u_{{1,1}}&u_{{1,2}}&u_{{1,3}}&\ldots&u_{{1,N_{y}}}\\ u_{{2,1}}&u_{{2,2}}&u_{{2,3}}&\ldots&u_{{2,N_{y}}}\\ \vdots\\ u_{{N_{x},1}}&u_{{N_{x},2}}&u_{{N_{x},3}}&\ldots&u_{{N_{x},N_{y}}}\end{pmatrix}.$

Zauważmy, że aby na przykład wyznaczyć iloraz różnicowy

$F=\dfrac{-u_{{i-1,j}}+2u_{{i,j}}-u_{{i+1,j}}}{h_{x}^{2}},$

nie musimy używać żadnej macierzy; możemy po prostu wykonać działanie, które, nieco nadużywając notacji Octave²²W Octave nie wolno indeksować macierzy od zera., zapiszemy

F = -U( 0:Nx-1, :) + 2*U( 1:Nx, :) - U( 2:Nx+1, :);
F = F/hx^2;

Analogicznie, dysponując tablicą Ax zawierającą wartości $a_{{i+1/2,j}}$ dla $i=0,\ldots,N_{x}$ oraz $j=1,\ldots,N_{y}$ , iloraz różnicowy

$F=\dfrac{-a_{{i-1/2,j}}u_{{i-1,j}}+(a_{{i-1/2,j}}+a_{{i+1/2,j}})u_{{i,j}}-a_{{i+1/2,j}}u_{{i+1,j}}}{h_{x}^{2}}$

wyznaczymy korzystając z operatorów tablicowego mnożenia,

F = -Ax(0:Nx-1).*U( 0:Nx-1, :) + (Ax(0:Nx-1)+Ax(1:Nx)).*U( 1:Nx, :) - Ax(1:Nx).*U( 2:Nx+1, :);
F = F/hx^2;

Powyższa obserwacja pozwoli nam efektywnie — bez użycia pętli i instrukcji warunkowych — obliczyć wynik działania dyskretnego dwuwymiarowego operatora dyfuzji na tablicy wartości funkcji siatkowej w węzłach.

Ćwiczenie 8.13

Zaimplementuj funkcję F = multdiff2dd(Uin, h, Ax, Ay, U), będącą analogonem opracowanej wcześniej funkcji jednowymiarowej. Niech $a(x,y)=xy+1$ . Wyznacz i zobrazuj wynik działania dyskretyzacji operatora dyfuzji $-\text{div}(a(x,y)\nabla u(x,y))$ w przypadku, gdy $u(x,y)=xy$ na $\bar{\Omega}=[0,2]\times[0,3]$ .

Rozwiązanie: [pokaż]

Najpierw implementacja operatora dyfuzji:

function F = multdiff2dd(Uin, h, A, U)
%
% wyznacza F = -div ( a(x,y) grad u(x,y) ), dla dyskretyzacji
%
% Uin - wartości u(x,y) we wnętrzu siatki
% h - długość kroku siatki: h = [hx,hy]
% A - wartości a(x,y) we wszystkich punktach pośrednich siatki: A = {Ax, Ay}
%	domyślnie a(x) = 1
% U - wartości zerowego rozszerzenia warunku brzegowego na wszystkie węzły siatki
%	domyślnie U(x) = 0

if nargin < 4 % domyślnie jednorodny warunek brzegowy
	U = zeros(size(Uin)+2);
end

[Nx, Ny] = size(U);

if nargin < 3 % domyślnie a(x) = 1
	Ax = ones(Nx-1,Ny-2);
	Ay = ones(Nx-2,Ny-1);
else
	Ax = A{1};
	Ay = A{2};
end

ix = 2:Nx-1; iy = 2:Ny-1; % indeksy odpowiadajace wartosciom z wnetrza
U(2:Nx-1,2:Ny-1) = Uin;
F = (-U(ix-1,iy).*Ax(ix-1,:) + U(ix,iy).*(Ax(ix-1,:)+Ax(ix,:)) - U(ix+1,iy).*Ax(ix,:))/(h(1)^2) + ...
	(-U(ix,iy-1).*Ay(:,iy-1) + U(ix,iy).*(Ay(:,iy-1)+Ay(:,iy)) - U(ix,iy+1).*Ay(:,iy))/(h(2)^2);
end

Podczas pisania tej funkcji musimy bardzo uważać, by uniknąć pomyłek przy przesunięciach indeksów. Zwróćmy także uwagę na pewną rozrzutność w gospodarowaniu pamięcią: wektory pomocnicze Ax i Ay są (formalnie niepotrzebnymi) kopiami pól struktury A: użyliśmy ich dla większej czytelności kodu.

Dalej, po podstawowych testach (które tutaj pominiemy wierząc, że Czytelnik samodzielnie i rutynowo je przeprowadzi), wykorzystamy świeżo napisaną funkcję w konkretnym przypadku opisanym w zadaniu.

X = 2; Y = 3;
Nx = 151; Ny = 153;
hx = X/(Nx+1); hy = Y/(Ny+1);

px = hx*[0:Nx+1]; py = hy*[0:Ny+1];

[pX, pY] = meshgrid(px,py);
A = pX.*pY+1;
contour(px, py, A ); title('A(x,y)'); pause(1);
A = permute(A, [2 1]); % przechodzimy na nasz układ tablicowy
Ax = 0.5*(A(1:end-1,2:end-1) + A(2:end,2:end-1));
Ay = 0.5*(A(2:end-1,1:end-1) + A(2:end-1,2:end));

U = permute(pX.*pY, [2 1]);
F = multdiff2dd(U(2:Nx+1,2:Ny+1),[hx,hy],{Ax,Ay},U);
contour(px(2:Nx+1), py(2:Ny+1), permute(F, [2 1]) ); title('-div (a(x,y) grad u(x,y))'); pause(1);

8.2.4. Rozwiązanie równania dyfuzji w obszarze dwuwymiarowym

Sama procedura rozwiązania równania nie zmieni się zanadto w stosunku do przypadku jednowymiarowego; dla ustalenia uwagi przyjmiemy, że $a(x,y)=xy+1$ oraz $f(x,y)=g(x,y)=1$ , a wszystkie pozostałe parametry zadania są takie same, jak w uprzednio rozważanym zagadnieniu Poissona w prostokącie.

%
% definicje siatki oraz macierzy Lx, Ly, Ix, Iy jak w przypadku r-nia Poissona -
% - tutaj pominięte!
%

global Lambda;
Lambda = kron(Iy,Lx) + kron(Ly,Ix);

% definicja wartości współczynnika dyfuzji w węzłach siatki
[pX, pY] = meshgrid(px,py);
A = pX.*pY+1;
A = permute(A, [2 1]);
Ax = 0.5*(A(1:end-1,2:end-1) + A(2:end,2:end-1));
Ay = 0.5*(A(2:end-1,1:end-1) + A(2:end-1,2:end));
% definicja operatorów niezbędnych do rozwiązania zadania
multdiff = @(x) reshape(multdiff2dd(reshape(x,Nx,Ny), [hx,hy], {Ax,Ay}), Nx*Ny, 1); % funkcja przykładająca operator dyfuzji
precond = @(x) reshape(fftsolve2dd(reshape(x, Nx, Ny)), Nx*Ny, 1); % jako operator ściskający wybieramy rozwiązanie r-nia Poissona metodą FFT

F = ones(Nx,Ny);
U = ones(Nx+2,Ny+2);
F = F - multdiff2dd(zeros(size(F)), [hx,hy], {Ax,Ay}, U); % uwzględnienie warunku brzegowego
[Uin, flag, relres, iter, resvec, eigest] = ...
		pcg(multdiff, F(:), (1e-3)*min([hx,hy])^2, 50, precond);
U(2:Nx+1,2:Ny+1) = reshape(Uin,Nx,Ny);

contour(px, py, permute(U, [2 1]) );
xlabel('x'); ylabel('y');

Brakującą funkcję fftsolve2dd bez trudu opracujesz, mając za wzór funkcję fftsolve3dd z przykładu 8.1.

Ćwiczenie 8.14 (trudne)

Zaimplementuj w Octave funkcję rozwiązującą trójwymiarowe zadanie dyfuzji.

Ćwiczenie 8.15 (trudne)

Zaimplementuj w Octave funkcję rozwiązującą na kostce $d$ -wymiarowe ( $d\in\{ 1,2,3\}$ ) zadanie dyfuzji z warunkiem brzegowym Neumanna.

8.3. Uwagi i uzupełnienia

W przypadku zadań trudniejszych, na przykład obszarów o bardziej skomplikowanej geometrii, warto skorzystać z bardziej wyspecjalizowanych narzędzi. Jedną z możliwości jest użycie biblioteki Deal.II napisanej w języku C++ i pozwalającej dyskretyzować równania różniczkowe różnych typów przy użyciu metody elementu skończonego. Zawiera on w sobie także kilka narzędzi do rozwiązywania zadań dyskretnych; ma także możliwość współpracy z innymi pakietami pomocniczymi, takimi jak zestaw procedur rozwiązywania zadań dyskretnych na komputerach wieloprocesorowych PETSc. Wadą tego podejścia jest konieczność nauczenia się i zrozumienia tej skomplikowanej biblioteki, co jest bardzo czasochłonne.

Tej wady nie mają proste narzędzia, takie jak program freeFEM++, ale ich stosowalność jest mocno ograniczona.

Ćwiczenie 8.16 (trudne)

Zaimplementuj w C/C++ funkcję rozwiązującą dwu- lub trójwymiarowe zadanie dyfuzji ze zmiennym współczynnikiem, z wykorzystaniem biblioteki Deal.II.

Ćwiczenie 8.17 (trudne)

Zaimplementuj w Octave lub w C++ funkcję rozwiązującą anizotropowe zadanie dyfuzji w $\Omega=(0,1)^{d}$ ze zmiennym współczynnikiem,

$-\text{div}(A(x)\nabla u(x))=f(x),\qquad x\in\Omega,$

z warunkiem brzegowym Dirichleta. Tutaj, dla zadanego $x\in\Omega$ , macierz $A(x)$ jest dodatnio określoną macierzą rozmiaru $d\times d$ . Rozważ $d\in\{ 1,2,3\}$ .

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

Obliczenia naukowe

Zagadnienia

8. Stacjonarne równanie dyfuzji

8.1. Laplasjan: dyskretyzacja metodą różnic skończonych

8.1.1. Dyskretyzacja jednowymiarowego laplasjanu

Ćwiczenie 8.1

8.1.2. Rozwiązanie równania Poissona w obszarze jednowymiarowym

Ćwiczenie 8.2

8.1.2.1. Testy

Ćwiczenie 8.3

8.1.3. Dyskretyzacja dwuwymiarowego laplasjanu

Ćwiczenie 8.4

8.1.4. Rozwiązanie równania Poissona w obszarze dwuwymiarowym

Ćwiczenie 8.5

8.1.5. Dyskretyzacja i rozwiązanie równania Poissona w obszarze trójwymiarowym

Ćwiczenie 8.6

8.1.5.1. Krytyka standardowego solvera bezpośredniego

8.1.5.2. Solver oparty na FFT

Przykład 8.1 (Trójwymiarowe równanie Poissona przez FFT)

Ćwiczenie 8.7

Quiz 8.1

Ćwiczenie 8.8

8.2. Równanie dyfuzji

8.2.1. Dyskretyzacja zadania jednowymiarowego

8.2.1.1. Uwzględnienie warunków brzegowych

Ćwiczenie 8.9

8.2.2. Rozwiązanie równania dyfuzji w obszarze jednowymiarowym

Ćwiczenie 8.10

Ćwiczenie 8.11

Ćwiczenie 8.12

8.2.3. Dyskretyzacja dwuwymiarowego operatora dyfuzji

Ćwiczenie 8.13

8.2.4. Rozwiązanie równania dyfuzji w obszarze dwuwymiarowym

Ćwiczenie 8.14 (trudne)

Ćwiczenie 8.15 (trudne)

8.3. Uwagi i uzupełnienia

Ćwiczenie 8.16 (trudne)

Ćwiczenie 8.17 (trudne)