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

©2024 mgr Jerzy Wałaszek
I LO w Tarnowie

Rozkład LU – macierz odwrotna

SPIS TREŚCI
Tematy pomocnicze

Problem

Dla danej macierzy An×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:

Xn×n =
 
x1,1 x1,2x1,n
x2,1 x2,2 x2,n
        
xn,1 xn,2 xn,n
 
   

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

X1 =
 
x1,1
x2,1
  
xn,1
 
   

Zgodnie z regułami mnożenia macierzy mamy:

A×
 
x1,1
x2,1

xn,1
 
   
=
 
1
0
…
0
 
   

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ń).

A×
 
x1,2
x2,2

xn,2
 
   
=
 
0
1
…
0
 
   

Postępując podobnie dla pozostałych kolumn X3, X4Xn, 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: Utwórz w X macierz jednostkową
K05: Dla i = 1,2,…,n: ; obliczamy kolejne kolumny
     wykonuj kroki K06…K07
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 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 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 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
Python (dodatek)
# Rozkład LU - macierz odwrotna
# Data: 28.04.2024
# (C)2024 mgr Jerzy Wałaszek
# -----------------------------

eps = 1e-12

# Funkcja dokonuje rozkładu LU macierzy A
#----------------------------------------
def ludist(A):
    n = len(A)
    for k in range(n-1):
        if abs(A[k][k]) < eps: return False
        for i in range(k+1,n):
            A[i][k] /= A[k][k]
        for i in range(k+1,n):
            for j in range(k+1,n):
                A[i][j] -= A[i][k]*A[k][j]
    return True

# Funkcja wyznacza wektor X
# na podstawie A i X[i]
#--------------------------
def lusolve(k,A,X):
    n = len(A)
    for i in range(1,n):
        s = 0.0
        for j in range(i):
            s += A[i][j]*X[j][k]
        X[i][k] -= s
    if abs(A[n-1][n-1]) < eps: return False
    X[n-1][k] /= A[n-1][n-1]
    for i in reversed(range(n-1)):
        s = 0.0
        for j in range(i+1,n):
            s += A[i][j]*X[j][k]
        if abs(A[i][i]) < eps: return False
        X[i][k] = (X[i][k]-s)/A[i][i]
    return True

# Program główny

# Odczytujemy wymiar macierzy A

n = int(input())

# Tworzymy macierze A i X

A = []
X = []

# Odczytujemy dane dla macierzy A

for i in range(n):
    a = input().split()
    a = [float(a[j]) for j in range(n)]
    A.append(a)
    X.append([0]*n)
    del a

# Wykonujemy rozkład LU macierzy A

if ludist(A):

    # Tworzymy macierz jednostkową w X   

    for i in range(n):
        X[i][i] = 1

# Wyznaczamy kolumny macierzy odwrotnej w X

    ok = True
    for i in range(n):
        if not lusolve(i,A,X):
            ok = False
            break
else: ok = False

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

if ok:
    for i in range(n):
        for j in range(n):
            print("%10.5f" % X[i][j],end="")
        print()
else: print("DZIELNIK ZERO")

# Usuwamy macierze z pamięci

del A,X
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
©2024 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.