disp 'metoda Choleskiego-Banachiewicza'

% rozklad LL' macierzy A (rozklad Choleskiego-Banachiewicza dla macierzy symetrycznej)

L = zeros(N); % utworzenie pustej macierzy L do wypelnienia

for i = 1 : N % dla kolejnych kolumn

% obliczenie elementow z diagonali

suma = 0;

for k = 1 : (i-1)

suma = suma + (L(i,k))^2;

end

L(i,i) = sqrt(A(i,i) - suma);

% obliczenie elementow spod diagonali

for j = (i+1) : N

suma = 0;

for k = 1 : (i-1)

suma = suma + L(j,k) * L(i,k);

end

L(j,i) = (A(j,i) - suma) / L(i,i);

end

end

L % wypisanie obliczonej macierzy L

LLT = L * L' % sprawdzenie czy L*L' == A

% wsteczne podstawienie (w dol) => wektor y = L' x

y = zeros(N,1); % utworzenie wektora y (chwilowo wypelnionego zerami)

for i = 1 : N

suma = 0;

for j = 1 : (i-1)

suma = suma + L(i,j) * y(j); % sumujemy iloczyny kolejnych y przez wspolczynniki

end

y(i) = ( b(i) - suma ) / L(i,i);

end

LT = L' % transpozycja L daje macierz gornotrojkatna

% ... dzieki temu mozna zastosowac wteczne podstawienie

% (w gore) => wektor x

x = zeros(N,1); % utworzenie wektora x

for i = N : -1 : 1

suma = 0;

for j = N : -1 : (i+1)

suma = suma + LT(i,j) * x(j);

end

x(i) = ( y(i) - suma ) / LT(i,i);

end

% wypisanie obliczonego x i sprawdzenie czy zostal dobrze policzony:

x_CB = x

Ax_b = A*x - b

%% a.2 wyznaczenie rozwiazania ukladu rownan metoda SOR

disp 'metoda SOR'

% parametr omega

w = 1

% rozklad macierzy A na L+D+U

L = tril(A,-1)

D = diag(diag(A))

U = triu(A,1)

% macierze B i c

B = inv(D + w*L) * ((1-w)*D - w*U)

c = w * inv(D + w*L) * b

X = []; % skasowanie zawartosci zmiennej X, jesli istniala

R = []; % skasowanie zawartosci zmiennej R, jesli istniala

X(:,1) = ones(N,1); % punkt startowy - przechowujemy jako pierwsza kolumne macierzy X

R(1) = norm(A*X(:,1)-b); % residuum dla punktu startowego - przechowujemy jako pierwszy element macierzy R

dokladnosc = 0.0001 % max norma wektora resztowego (residuum)

maxiter = 100 % max liczba iteracji

i = 1; % reset licznika iteracji

while (i <= maxiter) & (R(i) > dokladnosc)

X(:,i+1) = B * X(:,i) + c; % znalezienie i przechowanie w X kolejnego rozwiazania

R(i+1) = norm(A*X(:,i+1)-b); % znalezienie i przechowanie w R kolejnego residuum

i = i+1;

end

% wyswietlenie wynikow

x_SOR = X(:,i) % rozwiazanie koncowe

disp 'liczba iteracji ', i-1

[X' R'] % tabelka z wynikami i residuum

%% podpunkt b wyznacznik i macierz odwrotna do A

disp 'wyznacznik i macierz odwrotna do A'

detA = det(A)

Ainv = inv(A)