Serwis Edukacyjny
w I-LO w Tarnowie
obrazek

Materiały dla uczniów liceum

  Wyjście       Spis treści       Wstecz       Dalej  

Autor artykułu: mgr Jerzy Wałaszek

©2020 mgr Jerzy Wałaszek
I LO w Tarnowie

Rozkład LU – macierz odwrotna

SPIS TREŚCI
Tematy pomocnicze

Problem

Dla danej macierzy A n × n należy znaleźć macierz odwrotną wykorzystując rozkład LU macierzy A.

Rozwiązanie

Macierz odwrotna ( ang. the inverse of a matrix ) ( A -1 ) do macierzy A  spełnia równanie:

A  × A -1 = I  = A -1 × A
gdzie I  oznacza macierz jednostkową.

Oznaczmy X  = A -1. Wtedy powyższe równanie możemy zapisać jako:

A  × X = I

Niech:

Weźmy następnie pierwszą kolumnę macierzy X:

Zgodnie z regułami mnożenia macierzy mamy:

W wyniku otrzymujemy pierwszą kolumnę macierzy jednostkowej, gdyż z założenia X  = A -1. Wyrazy kolumny X 1 możemy potraktować jako niewiadome, a cały iloczyn jako układ n  równań liniowych względem tych niewiadomych. Układ ten rozwiązujemy metodą LU, którą opisaliśmy w poprzednim rozdziale. Następnie tak samo postępujemy z drugą kolumną X 2 ( rozkład LU wykonujemy tylko jeden raz na początku obliczeń – później wykorzystujemy go do rozwiązania kolejnych układów równań ).

Postępując podobnie dla pozostałych kolumn X 3, X 4 ... X n, wyznaczymy całą macierz X, która jest odwrotnością macierzy A. W algorytmie macierz X  może jednocześnie być macierzą jednostkową, której kolejne kolumny przekształcamy w kolumny X i. Macierz jednostkowa na końcu zostanie więc zastąpiona macierzą odwrotną do A.

Algorytm obliczania macierzy odwrotnej przy pomocy rozkładu LU

Wejście:

n  –  stopień macierzy, nN
A  – macierz kwadratowa stopnia n. Elementy A  ∈ R

Wyjście:

r = true, macierz kwadratowa stopnia n  X zawiera macierz odwrotną do A. Elementy X ∈ R
r  = false, błąd dzielenia przez zero

Zmienne pomocnicze:

i   –  indeks, iN

Lista kroków:

K01:  r  ← false zakładamy porażkę
K02: Wykonaj rozkład LU
dla macierzy A
 
K03: Jeśli rozkład się nie powiódł,
to zakończ
 
K04: Ustaw w X  macierz jednostkową  
K05: Dla i  = 1, 2, ...n:
wykonuj kroki K06...K07
obliczamy kolejne kolumny
K06:     Przyjmij za wyrazy wolne i-tą kolumnę macierzy X.
    Rozwiąż układ równań metodą LU.
    Zapisz wektor niewiadomych w i-tej kolumnie X
macierzy odwrotnej
K07:     Jeśli rozwiązanie się nie powiodło,
    to zakończ
 
K08:  r  ← true zaznaczamy sukces
K09: Zakończ  

Powyższy algorytm nie jest nowym algorytmem, wykorzystuje on dwa poprzednio opisane algorytmy rozkładu LU oraz rozwiązywania układu równań liniowych z wykorzystaniem macierzy L i U  otrzymanych w rozkładzie LU. Jedyną różnicą będzie odpowiednie potraktowanie kolumn macierzy X  jako wektora wyrazów wolnych B  oraz na wyjściu jako wektora niewiadomych. Zadanie to prosto zrealizujemy, modyfikując funkcję lusolve( ) tak, aby przyjmowała jako parametr numer kolumny macierzy X.


Przykładowe programy

Uwaga:

Zanim uruchomisz program, przeczytaj wstęp do tego artykułu, w którym wyjaśniamy funkcje tych programów oraz sposób korzystania z nich.

Program wczytuje dane definiujące macierz A  i wylicza macierz do niej odwrotną. Jeśli obliczenie nie może zostać wykonane, program wypisuje komunikat "DZIELNIK ZERO".

Dane dla macierzy A  są zdefiniowane następująco:

Pierwsza liczba określa stopień n.
Następne n  × n  liczb określa zawartość macierzy A  wierszami.

Dane przykładowe ( przekopiuj do schowka i wklej do okna konsoli ):

3
25 5 1
64 8 1
144 12 1
Pascal
// Rozkład LU - macierz odwrotna
// Data: 20.10.2010
// (C)2020 mgr Jerzy Wałaszek
//-----------------------------

program luinv;

type
  NReal = array of double; // typ tablic wierszy
  MReal = array of NReal;  // typ tablicy wskaźników

const
  eps = 1e-12;

// Funkcja dokonuje rozkładu LU macierzy A
//----------------------------------------
function ludist ( n : integer; A : MReal ) : boolean;
var
  i, j, k : integer;
begin
  for k := 0 to n - 2 do
  begin
    if abs ( A [ k ][ k ] ) < eps then Exit ( false );

    for i := k + 1 to n - 1 do
      A [ i ][ k ] := A [ i ][ k ] / A [ k ][ k ];

    for i := k + 1 to n - 1 do
      for j := k + 1 to n - 1 do
        A [ i ][ j ] := A [ i ][ j ] - A [ i ][ k ] * A [ k ][ j ];
  end;

  ludist := true;

end;

// Funkcja wyznacza wektor X na podstawie A i B
//---------------------------------------------
function lusolve ( k, n : integer; A, X : MReal ) : boolean;
var
  i, j : integer;
  s   : double;
begin
  for i := 1 to n - 1 do
  begin
    s := 0;

    for j := 0 to i - 1 do s := s + A [ i ][ j ] * X [ j ][ k ];

    X [ i ][ k ] := X [ i ][ k ] - s;
  end;

  if abs ( A [ n-1 ][ n-1 ] ) < eps then Exit ( false );

  X [ n-1 ][ k ] := X [ n-1 ][ k ] / A [ n-1 ][ n-1 ];

  for i := n - 2 downto 0 do
  begin
    s := 0;

    for j := i + 1 to n - 1 do s := s + A [ i ][ j ] * X [ j ][ k ];

    if abs ( A [ i ][ i ] ) < eps then Exit ( false );

    X [ i ][ k ] := ( X [ i ][ k ] - s ) / A [ i ][ i ];
  end;

  lusolve := true;

end;

// Program główny

var
  A, X  : MReal;
  n, i, j : integer;
  ok    : boolean;

begin

  // Odczytujemy wymiar macierzy A

  read ( n );

  // Tworzymy macierze A i X

  SetLength ( A, n );
  SetLength ( X, n );

  for i := 0 to n - 1 do
  begin
    SetLength ( A [ i ], n );
    SetLength ( X [ i ], n );
  end;   

  // Odczytujemy dane dla macierzy A

  for i := 0 to n - 1 do
    for j := 0 to n - 1  do read ( A [ i ][ j ] );

  // Wykonujemy rozkład LU macierzy A

  if ludist ( n, A ) then
  begin
  
  // Tworzymy macierz jednostkową w X    

    for i := 0 to n - 1 do
    begin
      for j := 0 to n - 1 do  X [ i ][ j ] := 0;
      X [ i ][ i ] := 1;
    end;

  // Wyznaczamy kolejne kolumny macierzy odwrotnej w X

    ok := true;
    for i := 0 to n - 1 do
      if not lusolve ( i, n, A, X ) then
      begin
        ok := false;
        break;
      end;
  end
  else ok := false;

  // Jeśli obliczenia się powiodły, wyświetlamy X

  if ok then
  begin
    for i := 0 to n - 1 do
    begin
      for j := 0 to n - 1 do write ( X [ i ][ j ] :10:5, ' ' );
      writeln;
    end;      
  end
  else writeln ( 'DZIELNIK ZERO' );

  // Usuwamy macierze z pamięci

  for i := 0 to n - 1 do
  begin
    SetLength ( A [ i ], 0 );
    SetLength ( X [ i ], 0 );
  end;
  SetLength ( A, 0 );
  SetLength ( X, 0 );
end.
C++
// Rozkład LU - macierz odwrotna
// Data: 20.10.2010
// (C)2020 mgr Jerzy Wałaszek
//-----------------------------

#include <iostream>
#include <iomanip>
#include <cmath>

using namespace std;

const double eps = 1e-12;

// Funkcja dokonuje rozkładu LU macierzy A
//----------------------------------------
bool ludist ( int n, double ** A )
{
  int i, j, k;

  for( k = 0; k < n - 1; k++ )
  {
    if( fabs ( A [ k ][ k ] ) < eps ) return false;

    for( i = k + 1; i < n; i++ )
      A [ i ][ k ] /= A [ k ][ k ];

    for( i = k + 1; i < n; i++ )
      for( j = k + 1; j < n; j++ )
        A [ i ][ j ] -= A [ i ][ k ] * A [ k ][ j ];
  }

  return true;
}

// Funkcja wyznacza wektor X na podstawie A i X [ i ] 
//---------------------------------------------------
bool lusolve ( int k, int n, double ** A, double ** X )
{
  int    i, j;
  double s;

  for( i = 1; i < n; i++ )
  {
    s = 0;

    for( j = 0; j < i; j++ ) s += A [ i ][ j ] * X [ j ][ k ];

    X [ i ][ k ] -= s;
  }

  if( fabs ( A [ n-1 ][ n-1 ] ) < eps ) return false;

  X [ n-1 ][ k ] /= A [ n-1 ][ n-1 ];

  for( i = n - 2; i >= 0; i-- )
  {
    s = 0;

    for( j = i + 1; j < n; j++ ) s += A [ i ][ j ] * X [ j ][ k ];

    if( fabs ( A [ i ][ i ] ) < eps ) return false;

    X [ i ][ k ] = ( X [ i ][ k ] - s ) / A [ i ][ i ];
  }

  return true;
}

// Program główny

int  main( )
{
  double **A, **X;
  int n, i, j;
  bool ok;

  cout << setprecision ( 5 ) << fixed;

  // Odczytujemy wymiar macierzy A

  cin >> n;

  // Tworzymy macierze A i X

  A = new double * [ n ];
  X = new double * [ n ];

  for( i = 0; i < n; i++ )
  {
    A [ i ] = new double [ n ];
    X [ i ] = new double [ n ];
  }

  // Odczytujemy dane dla macierzy A

  for( i = 0; i < n; i++ )
    for( j = 0; j < n; j++ ) cin >> A [ i ][ j ];

  // Wykonujemy rozkład LU macierzy A

  if( ludist ( n, A ) )
  {
  
  // Tworzymy macierz jednostkową w X    

    for( i = 0; i < n; i++ )
    {
      for( j = 0; j < n; j++ ) X [ i ][ j ] = 0;
      X [ i ][ i ] = 1;
    }

  // Wyznaczamy kolejne kolumny macierzy odwrotnej w X

    ok = true;
    for( i = 0; i < n; i++ )
      if( ! lusolve ( i, n, A, X ) )
      {
        ok = false;
        break;
      }
  }
  else ok = false;

  // Jeśli obliczenia się powiodły, wyświetlamy X

  if( ok )
  {
    for( i = 0; i < n; i++ )
    {
      for( j = 0; j < n; j++ )
        cout << setw ( 10 ) << X [ i ][ j ] << " ";
      cout << endl;
    }      
  }
  else cout << "DZIELNIK ZERO\n";

  // Usuwamy macierze z pamięci

  for( i = 0; i < n; i++ )
  {
    delete [ ] A [ i ];
    delete [ ] X [ i ];
  }
  delete [ ] A;
  delete [ ] X;

  return 0;
}
Basic
' Rozkład LU - macierz odwrotna
' Data: 20.10.2010
' (C)2020 mgr Jerzy Wałaszek
'-----------------------------

Const eps As Double = 1e-12

' Funkcja dokonuje rozkładu LU macierzy A
'----------------------------------------
Function ludist ( n As Integer, A As Double Ptr Ptr ) As Integer

  Dim As Integer i, j, k
  
  For k = 0 To n - 2
    If Abs ( A [ k ][ k ] ) < eps Then Return 0
    For i = k + 1 To n - 1: A [ i ][ k ] /= A [ k ][ k ] : Next
    For i = k + 1 To n - 1
      For j = k + 1 To n - 1: A [ i ][ j ] -= A [ i ][ k ] * A [ k ][ j ] : Next
    Next
  Next

  ludist = 1

End Function

' Funkcja wyznacza wektor X na podstawie A i X [ i ] 
'---------------------------------------------------
Function lusolve ( k As Integer, n As Integer, A As Double Ptr Ptr, X As Double Ptr Ptr ) As Integer

  Dim As Integer i, j 
  Dim As Double s

  For i = 1 To n - 1
    s = 0
    For j = 0 To i - 1: s += A [ i ][ j ] * X [ j ][ k ] : Next
    X [ i ][ k ] -= s
  Next

  If Abs ( A [ n - 1 ][ n - 1 ] ) < eps Then Return 0

  X [ n - 1 ][ k ] /= A [ n - 1 ][ n - 1 ] 

  For i = n - 2 To 0 Step -1
    s = 0
    For j = i + 1 To n - 1: s += A [ i ][ j ] * X [ j ][ k ] : Next
    If Abs ( A [ i ][ i ] ) < eps Then Return 0
    X [ i ][ k ] = ( X [ i ][ k ] - s ) / A [ i ][ i ] 
  Next

  lusolve = 1

End Function

Dim As Double Ptr Ptr A, X
Dim As Integer i, j, ok, n

Open Cons For Input As #1

' Odczytujemy wymiar macierzy A

Input #1, n

' Tworzymy macierze A i X

A = New Double Ptr [ n ] 
X = New Double Ptr [ n ] 
For i = 0 To n - 1
  A [ i ] = New Double [ n ] 
  X [ i ] = New Double [ n ] 
Next

' Odczytujemy dane dla macierzy A

For i = 0 To n - 1
  For j = 0 To n - 1
    Input #1, A [ i ][ j ] 
  Next
Next

Close #1

' Wykonujemy rozkład LU macierzy A

If ludist ( n, A ) = 1 Then
  
' Tworzymy macierz jednostkową w X    

  For i = 0 To n - 1
    For j = 0 To n - 1: X [ i ][ j ] = 0: Next
    X [ i ][ i ] = 1
  Next

' Wyznaczamy kolejne kolumny macierzy odwrotnej w X

  ok = 1
  For i = 0 To n - 1 
    If lusolve ( i, n, A, X ) = 0 Then
        ok = 0
        Exit For
    End If
  Next
Else
  ok = 0
End If

' Jeśli obliczenia się powiodły, wyświetlamy X

If ok = 1 Then
  For i = 0 To n - 1
    For j = 0 To n - 1: Print Using "####.##### ";X [ i ][ j ];: Next
    Print
  Next
Else
  Print "DZIELNIK ZERO"
End If 

' Usuwamy macierze dynamiczne

For i = 0 To n - 1
  Delete [ ] A [ i ] 
  Delete [ ] X [ i ] 
Next
Delete [ ] A
Delete [ ] X

End
Wynik:
3
25 5 1
64 8 1
144 12 1

   0.04762   -0.08333    0.03571
  -0.95238    1.41667   -0.46429
   4.57143   -5.00000    1.42857
Na początek:  podrozdziału   strony 

Zespół Przedmiotowy
Chemii-Fizyki-Informatyki

w I Liceum Ogólnokształcącym
im. Kazimierza Brodzińskiego
w Tarnowie
ul. Piłsudskiego 4
©2020 mgr Jerzy Wałaszek

Materiały tylko do użytku dydaktycznego. Ich kopiowanie i powielanie jest dozwolone
pod warunkiem podania źródła oraz niepobierania za to pieniędzy.

Pytania proszę przesyłać na adres email: i-lo@eduinf.waw.pl

Serwis wykorzystuje pliki cookies. Jeśli nie chcesz ich otrzymywać, zablokuj je w swojej przeglądarce.
Informacje dodatkowe.