-
Notifications
You must be signed in to change notification settings - Fork 0
/
rozkladCroutaUL.m
54 lines (54 loc) · 2.42 KB
/
rozkladCroutaUL.m
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
% (C) Jan Skwarek
function [U, L] = rozkladCroutaUL(A)
% sprawdzamy wymiary macierzy A
[rows, columns] = size(A);
% przypadek w ktorym macierz A nie jest kwadratowa - rozklad UL takiej
% macierzy nie istnieje - wyrzucamy blad
if rows ~= columns
error("Macierz nie jest kwadratowa!")
end
% standardowa implementacja algorytmu rozkladu Crouta, zmienione sa
% tylko nieco indeksy aby uzyskac rozklad UL zamiast rozkladu LU
for i = rows:-1:1
% uzupelniamy ostatni rzad macierzy L
L(rows, i) = A(rows, i);
% w rozkladzie Crouta macierz gornotrojkatna ma jedynki na glownej
% przekatnej
U(i, i) = 1;
end
for j = rows - 1:-1:1
if L(rows, rows) == 0 || abs(L(rows, rows)) < 1e-5
% wartosc w ostatnim wierszu w ostatniej kolumnie nie moze byc
% rowna zero - dzielenie przez zero - jeden z wiodacych minorow
% glownych ale patrzac od prawego dolnego rogu macierzy moze
% byc rowny zero
error("Dzielenie przez zero!" + ...
" Jeden z wiodacych minorow glownych, ale liczonych od" + ...
" prawego dolnego rogu macierzy wejsciowej jest rowny zero!" + ...
" Rozkład UL nie istnieje.")
end
% uzupelniamy ostatnia kolumne macierzy U
U(j, rows) = A(j, rows) / L(rows, rows);
end
for i = rows - 1:-1:1
for j = rows - 1:-1:i
% uzupelniamy kolejne kolumny macierzy L
L(j, i) = A(j, i) - dot(L(rows:-1:j + 1, i), U(j, rows:-1:j + 1));
end
for j = i - 1:-1:1
% sprawdzamy czy na glownej przekatnej macierzy L pojawilo sie
% zero. Jezeli tak, to przerywamy program i wyrzucamy blad -
% rozklad UL nie istnieje - jeden z wiodacych minorow glownych
% ale patrzac od prawego dolnego rogu macierzy wejsciowej A
% jest rowny zero
if L(i, i) == 0 || abs(L(i, i)) < 1e-5
error("det(L) bliski zeru. Dzielenie przez zero!" + ...
" Jeden z wiodacych minorow glownych, ale liczonych" + ...
" od prawego dolnego rogu macierzy wejsciowej jest rowny zero!" + ...
" Rozkład UL nie istnieje.")
end
% uzupelniamy kolejne kolumny macierzy U
U(j, i) = (A(j, i) - dot(L(rows:-1:i + 1, i), U(j, rows:-1:i + 1))) / L(i, i);
end
end
end