Rozkład LU – macierz odwrotna


Tematy pokrewne
Macierze
Podstawowe pojęcia dotyczące macierzy
Reprezentacja macierzy w pamięci komputera
Mnożenie macierzy przez skalar
Dodawanie macierzy
Transponowanie macierzy
Mnożenie macierzy
Potęgowanie macierzy
Eliminacja Gaussa
Eliminacja Gaussa-Crouta
Wyznacznik macierzy
Rozkład LU
Rozkład LU – rozwiązywanie układu równań liniowych
Rozkład LU – macierz odwrotna
Rozkład LU – wyznacznik macierzy
 

Problem

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

 

 

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 X1 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ą X2 (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 X3, X4 ... Xn 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 Xi. 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, n N
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.

Elementy pomocnicze:
i,  –  indeks i N
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 K06...K07 ; obliczamy kolejne kolumny
K06:     Przyjmij za wyrazy wolne kolumnę i-tą 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.

 

Program

Ważne:

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. Poniżej mamy przykładowe dane wejściowe dla programu:

 

Dane dla programu:

 

3
25 5 1
64 8 1
144 12 1

Lazarus
// Rozkład LU - macierz odwrotna
// Data: 20.10.2010
// (C)2012 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.
Code::Blocks
// Rozkład LU - macierz odwrotna
// Data: 20.10.2010
// (C)2012 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;
}
Free Basic
' Rozkład LU - macierz odwrotna
' Data: 20.10.2010
' (C)2012 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

 



List do administratora Serwisu Edukacyjnego Nauczycieli I LO

Twój email: (jeśli chcesz otrzymać odpowiedź)
Temat:
Uwaga: ← tutaj wpisz wyraz  ilo , inaczej list zostanie zignorowany

Poniżej wpisz swoje uwagi lub pytania dotyczące tego rozdziału (max. 2048 znaków).

Liczba znaków do wykorzystania: 2048

 

W związku z dużą liczbą listów do naszego serwisu edukacyjnego nie będziemy udzielać odpowiedzi na prośby rozwiązywania zadań, pisania programów zaliczeniowych, przesyłania materiałów czy też tłumaczenia zagadnień szeroko opisywanych w podręcznikach.



   I Liceum Ogólnokształcące   
im. Kazimierza Brodzińskiego
w Tarnowie

©2017 mgr Jerzy Wałaszek

Dokument ten rozpowszechniany jest zgodnie z zasadami licencji
GNU Free Documentation License.