Odwracanie macierzy
foster
Algorytm do odwracania macierzy przedstawiony poniżej pochodzi z książki Donalda E. Knuth'a "Sztuka programowania" tom 1. Odwraca on macierz bardzo efektywnie nie korzystając przy tym ze znanego z matematyki wzoru wykorzystującego wyznacznik macierzy i dopełnienia algebraiczne.
Algorytm zastepujący macierz jej macierzą odwrotną przy założeniu, że elementami macierzy są A[i,j] dla 1<=i,j<=n :
-
Dla k=1,2,...,n wykonuj co następuje: przejrzyj wszytkie kolumny wiersza k nie wykorzystane do tej pory jako kolumny osiowe w poszukiwaniu elementu o największej wartości bezwzglednej. Przypisz do C[k] numer kolumny, w której znaleziono taki element, i wykonaj operacje osiowania, przyjmując ten element za element osiowy. Jeśli wszystkie takie elementy są zerami, to macierz jest osobliwa i nie można jej odwrócić.)
-
Dokonaj permutacji wierszy i kolumn w ten sposob, żeby to, co było wierszem k, stało sie wierszem C[k], a to co było kolumną C[k], stało się kolumną k.
Operacja osiowania przebiega w następujący sposób:
Przed osiowaniem Po osiowaniu
kolumna inna kolumna inna
osiowa kolumna osiowa kolumna
/ . . \ / . . \
| . . | | . . |
| | | |
wiersz osiowy |.. a ......... b ..| |.. 1/a ...... b/a ...|
| | | |
| . . | | . . |
| . . | | . . |
| . . | | . . |
| | | |
|.. c ......... d ..| |..-c/a .... d-bc/a ..|
inny wiersz | | | |
| . . | | . . |
\ . . / \ . . /
Poniżej przedstawiam kod źródłowy programu realizujacego powyższy algorytm w języku Pascal.
{$N+}
Program Macierz_Odwrotna;
Const
Max = 50; { Maksymalny stopien macierzy }
Zero = 1.0E-7; { Umowne zero }
Type
TMacierz= Array[1..Max,1..Max] of Real;
Function OdwracanieMacierzy(Var Macierz: TMacierz;Stopien:Byte): Boolean;
Var
w,k,w1,k1,i, { zmienne wykorzystywane przy iteracji }
p,q : Integer; { wspolrzedne elementu osiowego }
c : Array[1..Max] of Byte;
BylaJuzKolOsiowa: Array[1..Max]of boolean;
Max : Real;
MacPom : TMacierz;
begin
FillChar(BylaJuzKolOsiowa,SizeOf(BylaJuzKolOsiowa),False);
for w:=1 to Stopien do
begin
c[w]:=0;
Max:=0.0;
for k:=1 to Stopien do
begin
if not(BylaJuzKolOsiowa[k]) then
if Max < Abs(Macierz[w,k]) then
begin
Max:=Abs(Macierz[w,k]);
p:=w;
q:=k;
end;
end;
if Max<Zero then
begin
OdwracanieMacierzy:= False;
Exit;
end;
for w1:=1 to Stopien do
for k1:=1 to Stopien do
if (w1<>p)and(k1<>q)then
Macierz[w1,k1]:= Macierz[w1,k1]-(Macierz[p,k1]*Macierz[w1,q]/Macierz[p,q]);
for w1:=1 to Stopien do
if w1<>p then Macierz[w1,q]:= -Macierz[w1,q]/Macierz[p,q];
for k1:=1 to Stopien do
if k1<>q then Macierz[p,k1]:= Macierz[p,k1]/Macierz[p,q];
Macierz[p,q]:= 1/Macierz[p,q];
c[w]:= q;
BylaJuzKolOsiowa[q]:= True;
end;
(* Permutowanie macierzy *)
MacPom:=Macierz;
for k:=1 to Stopien do
for i:=1 to Stopien do
Macierz[c[k],i]:= MacPom[k,i];
MacPom:=Macierz;
for k:=1 to Stopien do
for i:=1 to Stopien do
Macierz[i,k]:= MacPom[i,c[k]];
OdwracanieMacierzy:= True;
end;
Procedure Wczytaj(var Macierz:TMacierz;var Stopien:Byte);
Var
f : Text;
i,j : Integer;
Begin
Assign(f,'.\in.txt');
Reset(f);
ReadLn(f,Stopien);
if Stopien>max then
begin
WriteLn('Macierz za duza!!!');
Halt;
end;
for i:=1 to Stopien do
begin
for j:=1 to Stopien do
Read(f,Macierz[i,j]);
ReadLn(f);
end;
Close(f);
end;
Procedure Zapisz(Var Macierz:TMacierz; Stopien:Integer; OK:Boolean);
Var
f : Text;
i,j : Integer;
Begin
Assign(f,'.\out.txt');
Rewrite(f);
if OK then
for i:=1 to Stopien do
begin
for j:=1 to Stopien do
Write(f,Macierz[i,j],' ');
WriteLn(f);
end
else
WriteLn(f,'Macierz osobliwa!!!');
Close(f);
end;
Var
Macierz : TMacierz;
Stopien : Byte;
Begin
Wczytaj(Macierz,Stopien);
if OdwracanieMacierzy(Macierz,Stopien) then
Zapisz(Macierz,Stopien,True)
else
Zapisz(Macierz,Stopien,False);
end.
Program można ulepszyć pisząc dobra funkcję permutowania wierszy i kolumn, która nie będzie wykorzystywała macierzy pomocniczej. Stopień macierzy, którą można będzie odwrócić zwiększy się wtedy dwukrotnie (czyli do 100).
Metoda eliminacji Gaussa-Jordana Dane wejsciowe Dwie macierze -jedna odwracana,druga jednostkowa Na obu macierzach wykonujemy te same operacje Gdy pierwsza macierz sprowadzimy do macierzy jednostkowej druga przyjmie postac macierzy odwrotnej 1. Zamiana wierszy nie wplywa na macierz odwrotna 2. Dodanie wierszy nie wplywa na macierz odwrotna 3. Pomnozenie wierszy przez skalar rozny od zera nie wplywa na macierz odwrotna
uses
Crt;
{ Rozmiary macierzy. }
const
rozn = 20;
rozb = 40;
{ Macierz }
type
macierz = array [1..rozn,1..rozb] of real;
{ dane: }
{ A - macierz do odwr˘cenia (pierwszych n kolumn) }
{ n - rozmiar macierzy }
{ wynik: }
{ A - kolumny od n + 1 do 2n }
procedure odwmac (var a: macierz; n: integer);
var
i,j,k,m,wm: Integer;
s,t,w,p,max: Real;
b:macierz;
begin
for i:=1 to n do
for j:=1 to n do
b[i,j]:=a[i,j];
w:= 1;
for i:= 1 to n-1 do
begin
max:= abs (b[i,i]);
wm:= i;
for j:= i+1 to n do { wybor maxymalnego elementu }
if abs (b[j,i]) > max then
begin
max:= abs (b[j,i]);
wm:= j
end;
if wm > i then { zamiana wierszy }
begin
for k:= i to n do
begin
p:= b[i,k];
b[i,k]:= b[wm,k];
b[wm,k]:= p
end;
w:= -w { zamiana wierszy zmienia znak }
end;
if b[i,i] = 0 then { macierz osobliwa }
w:= 0;
if w <> 0 then { eliminacja Gaussa }
for j:= i+1 to n do
begin
p:= b[j,i] / b[i,i];
for k:= n downto i+1 do
b[j,k]:= b[j,k] - p * b[i,k]
end;
w:= w * b[i,i] { iloczyn elementow na przekatnej }
end;
w:= w * b[n,n];
if w=0 then
begin
writeln('Macierz osobliwa');
exit;
end;
for i := 1 to n do
for j := n + 1 to 2 * n do
if i = j - n then
A [i,j] := 1
else
A [i,j] := 0;
for i := 1 to n do
begin
max:= abs (a[i,i]);
wm:= i;
for j:= i+1 to n do { wybor maxymalnego elementu }
if abs (a[j,i]) > max then
begin
max:= abs (a[j,i]);
wm:= j
end;
if wm > i then { zamiana wierszy }
begin
for k:= i to 2*n do
begin
p:= a[i,k];
a[i,k]:= a[wm,k];
a[wm,k]:= p
end;
end;
s := A [i,i];
for j := 1 to 2 * n do
A [i,j] := A [i,j] / s;
for j := i + 1 to n do
begin
t := A [j,i];
for k := i to 2 * n do
A [j,k] := - A [i,k] * t + A [j,k];
end;
end;
for k := n + 1 to 2 * n do
begin
for i := n downto 1 do
begin
s := A [i,k];
for j := i + 1 to n do
s := s - A [i,j] * A [j,k];
A [i,k] := s;
end;
end;
end;
var
A: macierz;
s,t: Real;
j,i,n: Integer;
begin
ClrScr;
WriteLn ('Odwracanie macierzy');
WriteLn;
Write ('Podaj rozmiar macierzy n=');
ReadLn (n);
for i := 1 to n do
begin
WriteLn ('Wprowad« ',i,'-ty wiersz macierzy');
for j := 1 to n do
ReadLn (A [i,j]);
end;
odwmac (A,n);
for i := 1 to n do
begin
for j := n + 1 to 2 * n do
Write (A [i,j]:4:8,' ');
WriteLn;
end;
ReadKey;
end.
U Bernarda Barona (Metody Numeryczne w Pascalu , kody źródłowe powinny być jeszcze dostępne)
można znaleźć procedurę odwracającą macierz bez dołączania macierzy jednostkowej
Uszanowanie dla autora. Algorytm oczywiście działa. Przeliczyłem nie macierz a potem sprawdziłem dwoma sposobami. Mam prośbę, moglibyście wyjaśnić dlaczego to działa?
Próbowałem z książki Knuth-a ale nie zrozumiałem jaki jest związek z "przeosiowaniem macierzy" a jej macierzą odwrotną.Pozdrawiam serdecznie
Nawet fajny algorytm. Tylko ciekawy jestem czy jest on zawsze stabilny. Metoda Gaussa-Jordana moze nie jest najszybsza ale zawsze doprowadzi nas do poprawnego wyniku. Oczywiscie gdy macierz jest nie osobliwa :)
Ładnie, nie czytałem tekstu tylko kod więc jeśli to Twój algorytm to bardzo ładny. Chłopaki - bez przesady Guass był matematykiem a foster jest nie ujmując mu ambitnym programistą, więc o co chodzi ;)
Gauss na pewno jest prostszy, a co za tym idzie, w praktyce moze okazac sie szybszy, zwlaszcza dla mniejszych macierzy (stale pomijane przy szacowaniu zlozonosci, czesto nie sa takie male jakby sie wydawac moglo). A jeszcze lepiej Gaussa-Jordana...
Ale swoja droga, to ciekawa metoda :)
A czy przypadkiem algorytm Gaussa nie jest szybszy?