Dzieląc równanie przez stałą dostajemy (8.20) dla prawej strony równej . Ponieważ załamuje się w , weźmiemy siatkę, w której jeden z węzłów wypadnie akurat w . W przypadku siatki jednorodnej zagwarantujemy to sobie, biorąc nieparzyste . Wstępnie wydaje się, że biorąc , uzyskamy dostatecznie dużą rozdzielczość schematu (wszak wówczas ). 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 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 (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 , 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 . 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,

Zgodnie z wyżej opisanym schematem, musimy wyznaczyć macierz odpowiadającą równaniu różnicowemu

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

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 — nie będzie klasyczne. Owszem, możemy spodziewać się rozwiązań słabych — klasy — 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 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 ? Przyczyną jest struktura, jaką mają macierze pX, pY, które Octave generuje poleceniem meshgrid. Aby ułatwić wizualizację takich macierzy, indeksy odpowiadające zmiennej 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.

Tak sformułowane zadanie ma charakter czysto jakościowy, a nie ilościowy — rzeczywiście, skoro 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 ! 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 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.

Faktycznie, mamy

A więc można (i nawet: trzeba!) w tym wypadku użyć solvera bezpośredniego. My implementujemy solver iteracyjny jedynie w celu oswojenia i sprawdzenia podejścia związanego z wykorzystaniem metody iteracyjnej…

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

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