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

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,

gdzie ui ma odpowiadać wartości rozwiązania w węźle dyskretyzacji xi, ui≈u⁢xi, oraz xi=i⋅hx, przy czym hx=X/Nx+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 UNx=u1,…,uNxT,

gdzie

Ponieważ interesują nas dobre przybliżenia (to znaczy przypadek, gdy hx jest małe — lub bardzo małe), znaczy to, że Nx będzie (bardzo) duże. Nasza macierz jest macierzą rozrzedzoną, bo ma nieco mniej niż 3⁢Nx niezerowych elementów. Jej współczynnik wypełnienia (density), czyli stosunek liczby niezerowych do liczby wszystkich elementów macierzy, rzędu 1/Nx, maleje ze wzrostem Nx. Nie ma więc sensu pamiętanie wszystkich Nx2 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);

Dzięki reprezentacji naszej macierzy w formacie rzadkim, macierz rozmiaru Nx zajmie nam w pamięci tylko około 8⋅3⋅Nx=24⁢⁢Nx 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⁢⁢Nx2 bajtów w przypadku, gdybyśmy reprezentowali ją jako macierz pełną.

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 L⁢D⁢LT) pozwoli nam wyznaczyć rozwiązanie układu TNx⁢UNx=FNx kosztem O⁢Nx 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;

Chcąc rozwiązać dwuwymiarowe zadanie Poissona,

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

Będziemy zatem mieli do czynienia z wartościami rozwiązania w punktach xi,yj, gdzie

przy czym

Wartości rozwiązania dyskretnego ui,j będą odpowiadać wartościom rozwiązania w węłach dyskretyzacji xi,yj, tzn. ui,j≈u⁢xi,yj.

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:

i analogicznie dla ∂2∂⁡y2⁢xi,yj. Otrzymamy więc w końcu równanie różnicowe postaci

gdzie fi,j=f⁢xi,yj. Mamy więc do wyznaczenia N=Nx⋅Ny niewiadomych, które w naturalny sposób układają się w macierz

(zwróćmy uwagę na to, że ,,fizyczne” współrzędne są odwrócone: wartości odpowiadające kolejnym współrzędnym w kierunku osi O⁢X 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.:

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

gdzie PN jest pięciodiagonalną macierzą o blokowej strukturze,

W kolejnych diagonalnych blokach PN znajdują się macierze trójdiagonalne TNx+2hy2⁢I, gdzie TNx jest znaną już nam macierzą dyskretyzacji jednowymiarowego laplasjanu. Jest tych bloków Ny. Używając symbolu Kroneckera dla iloczynu tensorowego macierzy mamy po prostu

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

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 PN 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);

W przypadku macierzy PN, prosta strategia polegająca na eliminacji Gaussa (z wykorzystaniem dodatkowo faktu, że PN jest pasmowa, o szerokości pasma 2⁢Nx), dałaby nam szansę rozwiązać to równanie kosztem O⁢Nx2⁢⁢Ny2, co jest, przy dużych wielkościach Nx i Ny, 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;

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

Czytelnik zechce sam sprawdzić, że dyskretyzacja operatora Laplace'a w kostce Ω=0,X×0,Y×0,Z różnicami skończonymi na regularnej siatce o oczku hx,hy,hz i Nx×Ny×Nz wewnętrznych węzłach, prowadzi do zadania różnicowego

dla i=1,…,Nx,⁢j=1,…,Ny,⁢k=1,…,Nz, z odpowiednimi warunkami brzegowymi. Oczywiście ui,j,k≈u⁢i⁢hx,j⁢hy,k⁢hz. Mamy więc do wyznaczenia N=Nx⋅Ny⋅Nz niewiadomych wewnątrz obszaru, które w naturalny sposób układają się w trójwymiarową tablicę

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:

Będziemy szukać wektora UN, który spełnia liniowy układ równań

gdzie SN jest siedmiodiagonalną macierzą o blokowej strukturze,

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 Nx=Ny=Nz=100 dochodzimy do N=106 (!) Dodatkowo, struktura naszej (bardzo rozrzedzonej) macierzy jest na tyle skomplikowana dla bezpośredniego solvera opartego na eliminacji Gaussa, że szybko staje się on nieskuteczny.

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

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

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

gdzie za ai+1/2 można wziąć na przykład

gdzie, zgodnie z intuicją, xi+1/2=xi+xi+12. Widzimy więc, że obliczenie działania dyskretnego operatora dyfuzji DNx:RNx+2→RNx na wektorze UNx (dane wzorem (8.20)), wymaga znajomości wartości także UNx na brzegu obszaru i prowadzi do wyznaczenia wartości wyniku jedynie w węzłach wewnętrznych.

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;

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

z warunkiem brzegowym

Mamy do wyznaczenia N=Nx⋅Ny niewiadomych, które w naturalny sposób układają się w macierz

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

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;

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.

Sama procedura rozwiązania równania nie zmieni się zanadto w stosunku do przypadku jednowymiarowego; dla ustalenia uwagi przyjmiemy, że a⁢x,y=x⁢y+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.