Algorytmy

Odwracanie macierzy

  • 2007-09-09 02:08
  • 8 komentarzy
  • 4932 odsłony
  • Oceń ten tekst jako pierwszy
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 :

1)  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ć.)
2) 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.

8 komentarzy

Brak avatara
Maciej_ZB 2014-01-29 13:30

Wyniki obliczania procedurą odwmac nie są takie same jak wyniki przeliczania macierzy na stronie http://calcoolator.pl/kalkulator_macierzy.html. Ktoś źle liczy

Brak avatara
Maciej_ZB 2014-01-29 11:51

Procedura odwmac ma błąd: dzielenie przez zero. Przy ostatniej (n-tej) iteracji w pętli for zmienna s przyjmuje wartość 0 (linia kodu: s := A [i,i];), a dwie linie dalej jest pętla for z dzieleniem przez s (linia kodu: A [i,j] := A [i,j] / s;).

mjankow 2009-11-07 20:26

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

bociang 2004-01-31 23:58

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

Rabbitsoft 2004-01-22 11:57

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

Dryobates 2004-01-20 23:45

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

Marooned 2004-01-20 22:15

A czy przypadkiem algorytm Gaussa nie jest szybszy?

my_nick 2004-01-20 19:27

1. Tag < delphi >
2. Formatowanie kodu (wszystkie słowa kluczowe - function, procedure, var, begin - małą literą i wcięcia)
3. Ten "rysunek" na początku się ciut rozpadł