Algorytmy i Struktury Danych - Rozkład LU

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:

X_n×n =

	x_1,1x_1,2 … x_1,nx_2,1x_2,2 … x_2,n ……… … x_n,1x_n,2 … x_n,n

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

X₁ =

	x_1,1x_2,1 … x_n,1
	x_1,1x_2,1 … x_n,1

Zgodnie z regułami mnożenia macierzy mamy:

A×

	x_1,1x_2,1 … x_n,1
	x_1,1x_2,1 … x_n,1

	1 0 …0
	1 0 …0

W wyniku otrzymujemy pierwszą kolumnę macierzy jednostkowej, gdyż z założenia X = A^-1. Wyrazy kolumny X₁ 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₂ (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×

	x_1,2x_2,2 … x_n,2
	x_1,2x_2,2 … x_n,2

	0 1 …0
	0 1 …0

Postępując podobnie dla pozostałych kolumn X₃, X₄ … 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; 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):

325 5 1
64 8 1144 12 1

Pascal

// Rozkład LU-macierz odwrotna
// Data: 20.10.2010
// (C)2020 mgr Jerzy Wałaszek
//---------------------------

program luinv;

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

// Przybliżenie zera
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):
    s = input().split()
    arr = [float(s[j]) for j in range(n)]
    a.append(arr)
    x.append([0]*n)
# 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")

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

do podrozdziału do strony

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.

Rozkład LU – macierz odwrotna