ROZWI
ZYWANIE UKA
ADU
ROZWI
ZYWANIE UKA
ADU
RÓWNAC
LINIOWYCH
RÓWNAC
LINIOWYCH
METOD
ROZKA
ADU
METOD
ROZKA
ADU
NA MACIERZE TRÓJK
TNE
NA MACIERZE TRÓJK
TNE
J. Chudzicka 1
Opis metody LU
Metoda LU jest metodą rozwiązywania układu
równań liniowych postaci:
W zapisie macierzowym A*x = y, gdzie A - macierz
współczynników, x - wektor niewiadomych, y - wektor
wyrazów wolnych.
J. Chudzicka 2
Pozwala także na szybkie wyliczenie wyznacznika
macierzy A .
W metodzie LU macierz współczynników
zapisywana jest jako iloczyn dwóch macierzy
trójkątnych: dolnej (ang. lower - L) i górnej (ang.
upper - U), tj. z elementami zerowymi - odpowiednio
- powyżej i poniżej przekątnej macierzy:
A = L*U
J. Chudzicka 3
Układ równań przyjmuje wówczas postać: L*U*x = y,
a jego rozwiÄ…zanie sprowadza siÄ™ do rozwiÄ…zania
dwóch układów równań z macierzami trójkątnymi,
które z kolei rozwiązuje się bardzo prosto:
L*z = y,
U*x = z
J. Chudzicka 4
Ostatecznie liczba mnożeń, potrzebnych do
wyznaczenia wektora x, wynosi n2, a dodawań n2 - n.
Zalety metody:
Øð
bardzo oszczędna gospodarka pamięcią;
Øð
wymaga najmniejszej liczby operacji w
porównaniu z innymi metodami dokładnymi (nie
biorÄ…c pod uwagÄ™ procedur specjalnych).
J. Chudzicka 5
Rozkład LU
Podstawowym problemem numerycznym w tej metodzie
jest dokonanie rozkładu LU macierzy współczynników.
Żeby ten rozkład macierzy był jednoznaczny zakłada się,
że elementy na głównej przekątnej jednej z macierzy, L
albo U, są równe 1.
Rozkład LU jest wyznaczany za pomocą metody
Doolittle'a.
Ta metoda nie jest niezawodna, tzn. podczas obliczeń
może wystąpić dzielenie przez zero. Istnieje jej
modyfikacja pozbawiona tej wady, nazywana metodÄ…
Doolittle-Crouta, w której wykorzystuje się częściowy
wybór elementu podstawowego.
J. Chudzicka 6
Element podstawowy to taki element w macierzy A, który
jest używany do rugowania zmiennych (czyli zerowania
odpowiadających im współczynników) z kolejnych równań.
Przy stosowaniu metody Doolittle'a wybiera siÄ™ element
podstawowy zawsze z przekątnej głównej i jeśli akk jest
równe zero, metoda zawodzi.
W metodach zmodyfikowanych wybierany jest ten element z
danej k-tej kolumny, który ma największy moduł. Następnie
wiersz, w którym znajduje się wybrany element, zamieniany jest
z k-tym wierszem, co powoduje, że element podstawowy
pojawia się na przekątnej głównej.
To gwarantuje, że podczas obliczeń nie wystąpi dzielenie przez
zero.
J. Chudzicka 7
Metoda Doolittle'a
W metodzie tej równość A = L*U traktuje
się jako układ n2 równań z n2
niewiadomymi. Te niewiadome to elementy
lij dla i > j (elementy poniżej przekątnej),
oraz uij dla j >= i (elementy na i powyżej
przekątnej), przy założeniu, że na diagonali
macierzy L znajdujÄ… siÄ™ 1:
J. Chudzicka 8
Wyznaczanie kolejnych elementów macierzy L i U
robi siÄ™ naprzemiennie, tj. raz wyznacza wiersz
macierzy U, raz kolumnÄ™ macierzy L.
J. Chudzicka 9
Wzory ogólne na poszczególne elementy macierzy
rozkładu przedstawiają się następująco:
Dla wszystkich
dla
dla
Z ostatniego równania wynika, że metoda
nie zadziała, gdy uii = 0.
J. Chudzicka 10
Liczba działań potrzebna do rozkładu:
mnożenia:
dodawania:
J. Chudzicka 11
Przykład (macierz 3x3)
J. Chudzicka 12
Pierwszy wiersz macierzy U:
J. Chudzicka 13
Pierwsza kolumna macierzy L:
J. Chudzicka 14
Drugi wiersz macierzy U:
J. Chudzicka 15
Druga kolumna macierzy L:
J. Chudzicka 16
Trzeci wiersz macierzy U:
J. Chudzicka 17
Metoda Banachiewicza
W metodzie Banachiewicza rozwiÄ…zywania
układów równań liniowych jest wykorzystywany
rozkład macierzy współczynników A na iloczyn
macierzy trójkątnej dolnej L i macierzy
trójkątnej górnej U  mającej tę własność, że
wszystkie elementy jej przekątnej głównej są
jedynkami (poprzednio dolna macierz zawierała
1 na przekÄ…tnej).
A = L*U
J. Chudzicka 18
FRAGMENTY KODU W DELPHI 7
Opis obiektów:
unit Obliczenia; &
const nmax = 20;
type
Tabl = array[1..nmax,1..nmax] of Real;
var
Form3: TForm3;
i,j,k,n,m: Integer;
nazwa: String;
A,B,L,U: Tabl;
Wczytywanie da
plik: Text;
det: Real;
J. Chudzicka 19
procedure Rozklad_LU(n: Integer; A: Tabl; var L,U: Tabl);
var
i,j,k: Integer;
sp: Real;
begin
for i:=1 to n do
for j:=1 to n do begin
if j=1 then L[i,j]:=A[i,j]
else L[i,j]:=0;
if j=i then U[i,j]:=1
else U[i,j]:=0;
end;
if Abs(L[1,1])<1e-10 then Halt(1);
for j:=2 to n do
Jeśli na przekątnej
U[1,j]:=A[1,j]/L[1,1];
b. małe wartości, to
for i:=2 to n do
zatrzymanie
L[i,2]:=A[i,2]-L[i,1]*U[1,2];
obliczeń
if Abs(L[2,2])<1e-10 then Halt(1);
J. Chudzicka 20
for j:=3 to n do begin
for i:=2 to j -1 do begin
sp:=0;
for k:=1 to i - 1 do
sp:=sp+L[i,k]*U[k,j];
U[i,j]:=(A[i,j]-sp)/L[i,i];
end;
for i:=j to n do begin
sp:=0;
for k:=1 to j-1 do
sp:=sp+L[i,k]*U[k,j];
L[i,j]:=A[i,j]-sp;
if (i=j) and (Abs(L[i,i])<1e-10) then Halt(1);
end;
end;
end;
Koniec procedury Rozklad_LU
J. Chudzicka 21
procedure MetBanach(n,m: Integer; L,U: Tabl; var B: Tabl;
var det: Real);
var
i,j,k: Integer;
sp: Real;
Y: Tabl;
begin
for k:=1 to m do
RozwiÄ…zanie
Y[1,k]:=B[1,k]/L[1,1];
pierwszego
for i:=2 to n do
układu
for k:=1 to m do begin
trójkątnego
sp:=0;
wyznaczajÄ…cego
for j:=1 to i-1 do
wektor Y
sp:=sp+L[i,j]*Y[j,k];
Y[i,k]:=(B[i,k]-sp)/L[i,i];
end;
J. Chudzicka 22
 for k:=1 to m do
RozwiÄ…zanie
B[n,k]:=Y[n,k];
drugiego układu
for i:= n -1 downto 1 do
trójkątnego
for k:=1 to m do begin
wyznaczajÄ…cego
sp:=0;
wektor X i
for j:=i+1 to n do
zapamiętanie go
sp:=sp+U[i,j]*B[j,k];
w macierzy B
B[i,k]:=Y[i,k]-sp;
end;
det:=1;
for i:=1 to n do
Obliczanie wyznacznika
det:=det*L[i,i];
end;
Koniec procedury MetBanach
J. Chudzicka 23
procedure TForm3.FormCreate(Sender: TObject);
begin
nazwa:='Pr_2_2.dan';
Dane.Lines.Text:='';
end;
procedure TForm3.BitBtn1Click(Sender: TObject);
begin
Form2.Show;
AssignFile(plik,nazwa);
Reset(plik); Readln(plik,n,m);
for i:=1 to n do begin
Wczytywanie
for j:=1 to n do
wierszami wraz z
Read(plik,A[i,j]);
macierzÄ…
for k:=1 to m do
wyrazów wolnych
Read(plik,B[i,k]);
end;
CloseFile(plik);
J. Chudzicka 24
AssignFile(plik,Edit1.Text);
Rewrite(plik); Writeln(plik,'PROGRAM 2.2.');
Writeln(plik,'Rozwiązywanie układu równań liniowych.');
Writeln(plik,'Metoda Banachiewicza.'); Writeln(plik);
Writeln(plik,'Liczba równań układu - n =',n:3);
Writeln(plik,'Liczba prawych stron - m =',m:3); Writeln(plik);
Writeln(plik,'Macierz współczynników:');
for i:=1 to n do begin
Writeln(plik,' wiersz nr ',i:3);
Wydruk
k:=0; Write(plik,' ');
macierzy A
for j:=1 to n do begin
po 5
k:=k+1;
elementów
if k=5 then begin
w wierszu
k:=0; Writeln(plik);
Write(plik,' ');
end;
Write(plik,' ',A[i,j]:16);
end; Writeln(plik);
end;
J. Chudzicka 25
Writeln(plik);
Writeln(plik,'Wektory prawych stron:');
for i:=1 to n do begin
Writeln(plik,' wiersz nr ',i:3);
Write(plik,' ');
for k:=1 to m do
Write(plik,' ',B[i,k]:16);
Writeln(plik);
end;
Writeln(plik);
Wydruk macierzy wyrazów wolnych (wierszami)
J. Chudzicka 26
Rozklad_LU(n,A,L,U); Wywołanie proc. Rozklad_LU
Wywołanie proc. MetBanach
MetBanach(n,m,L,U,B,det);
Writeln(plik,'Wyznacznik - det = ',det:16);
Writeln(plik);
Writeln(plik,'Rozwiązania układów równań:');
for i:=1 to n do begin
Writeln(plik,' wiersz nr ',i:3);
Write(plik,' ');
for k:=1 to m do
Write(plik,' ',B[i,k]:16);
Writeln(plik);
end;
Writeln(plik);
Za pomocÄ… procedury Rozklad_LU wyznaczane sÄ… macierze
L i U
J. Chudzicka 27
Wyznaczanie macierzy odwrotnej
Z definicji macierzy odwrotnej X wiadomo, że dla macierzy
kwadratowej nieosobliwej A (tzn. istnieje dla niej macierz
odwrotna) mamy zależność:
a11 a12 & a1n x11 x12 & x1n
1 0 & 0
a21 a22 & a2n x21 x22 & x2n 0 1 & 0
=
& & & & ..
& & & & & & & & & &
0 0 & 1
an1 an2 & ann xn1 xn2 & xnn
J. Chudzicka 28
Mnożąc macierze A i A-1 uzyskujemy n układów
równań względem n2 niewiadomych xij
a11 a12 & a1n x1j 0
0
a21 a22 & a2n x2j
&
& & & & & .. &
j = 1,2, & , n
=
1
aj1 aj2 & ajn Xjj
Wartość 1 jest
&
tylko w wierszu
& & & & & . &
o numerze j
0
an1 an2 & ann
xnj
J. Chudzicka 29
Tworzenie macierzy wyrazów
wolnych B dla n układów
for i:=1 to n do
for j:=1 to n do
if i=j then B[i,j]:=1
else B[i,j]:=0;
MetBanach(n,n,L,U,B,det);
Rozkład na macierze trójkątne już jest, więc od
razu wywołano procedurę MetBanach
J. Chudzicka 30
 Writeln(plik,'Macierz odwrotna:');
for i:=1 to n do begin
Writeln(plik,' wiersz nr ',i:3);
k:=0; Write(plik,' ');
for j:=1 to n do begin
k:=k+1;
if k=5 then begin
k:=0; Writeln(plik);
Write(plik,' ');
end;
Write(plik,' ',B[i,j]:16);
end;
Writeln(plik);
end;
Writeln(plik); CloseFile(plik);
Form2.Wyniki.Lines.LoadFromFile(Edit1.Text);
end;
Koniec procedury TForm3.BitBtn1Click ze slajdu 24
J. Chudzicka 31
procedure TForm3.BitBtn2Click(Sender: TObject);
begin
Form1.Close;
Form2.Close;
Form3.Close;
Form4.Close;
Form5.Close;
end;
procedure TForm3.BitBtn3Click(Sender: TObject);
begin
if OpenDialog1.Execute then begin
nazwa:=OpenDialog1.FileName;
Dane.Lines.LoadFromFile(nazwa);
end;
end;
J. Chudzicka 32
procedure TForm3.BitBtn4Click(Sender: TObject);
begin
Close;
end;
procedure TForm3.SpeedButton1Click(Sender: TObject);
begin
if SaveDialog1.Execute then begin
Dane.Lines.SaveToFile(SaveDialog1.FileName);
nazwa:=SaveDialog1.FileName;
end;
end;
end. Koniec programu
J. Chudzicka 33