Algorytmy i Struktury Danych - Rozkład LU-wyznacznik

Pamiętamy, iż algorytm dotyczący obliczania wyznacznika metodą rozwinięcia Laplace'a, miał klasą złożoności obliczeniowej O(n!), czyli bardzo niekorzystną. Na szczęście istnieją dużo szybsze algorytmy, realizujące to zadanie. Po pierwsze zauważmy, iż wyznacznik macierzy trójkątnej jest równy iloczynowi elementów leżących na głównej przekątnej. Można to bardzo prosto wyprowadzić z rozwinięcia Laplace'a. Pomiędzy macierzami A, L i U zachodzą odpowiednie zależności:

A = L×U
det A = det(L×U)
det A = det L×det U

det L = 1, ponieważ jest macierzą trójkątną, 
          a wszystkie elementy przekątnej są równe 1 
det A = 1×det U = det U

Zatem w celu obliczenia wyznacznika macierzy A wystarczy dokonać jej rozkładu LU i obliczyć iloczyn elementów leżących na przekątnej macierzy U.

Algorytm obliczania wyznacznika za pomocą rozkładu LU

Wejście:

n : stopień macierzy; n ∈ N.
A : macierz kwadratowa stopnia n; elementy A ∈ R.

Wyjście:

Wartość wyznacznika w zmiennej det lub informacja o błędzie; det ∈ R.

Elementy pomocnicze:

i : indeks; i ∈ N.

Lista kroków:

K01: Wyznacz rozkład LU macierzy A, 
     wynik w A
K02: Jeśli się powiodło, 
     to idź do kroku K05
K03: Pisz "DZIELNIK ZERO"
K04: Zakończ
K05: det ← A[1, 1] ; liczymy wyznacznik
K06: Dla i = 2, 3, …, n: ; kolejne iloczyny wyrazów
     wykonuj det ← det×A[i, i] ; na przekątnej głównej
K07: Pisz det
K08: Zakończ

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 odczytuje ze standardowego wejścia macierz A (najpierw liczba n określająca stopień macierzy, a następnie n² liczb dla kolejnych elementów), wyznacza jej rozkład LU za pomocą algorytmu Doolitle'a i jeśli się powiódł, oblicza i wyświetla wyznacznik macierzy A.

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

3
1 2 3
6 5 4
3 7 2

Pascal

// Wyznacznik przez rozkład LU
// Data: 10.02.2011
// (C)2020 mgr Jerzy Wałaszek
//-----------------------------

program ludet;

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

const
  eps = 1e-12;

// Funkcja realizuje algorytm
// rozkładu LU
//---------------------------
function lu(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;
  lu := true;
end;

// Program główny
//---------------

var
  A     : MReal;
  n, i, j : integer;
  det   : double;

begin

  // odczytujemy stopień macierzy A
  read(n);

  // tworzymy macierz A o odpowiednich rozmiarach
  SetLength(A, n);
  for i := 0 to n-1 do
    SetLength(A[i], n);

  // odczytujemy dane dla macierzy A
  for i := 0 to n-1 do
    for j := 0 to n-1 do
      read(A[i][j]);

  // obliczamy wyznacznik
  if lu(n, A) then
  begin
    det := A[0][0];
    for i := 1 to n-1 do
      det := det*A[i][i];
    writeln(det:0:4);
  end
  else
    writeln('DZIELNIK ZERO');

  // usuwamy macierz z pamięci
  for i := 0 to n-1 do
    SetLength(A[i], 0);
  SetLength(A, 0);

end.

C++

// Wyznacznik przez rozkład LU
// Data: 10.02.2011
// (C)2020 mgr Jerzy Wałaszek
//-----------------------------

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

using namespace std;

const double eps = 1e-12;

// Funkcja realizuje algorytm rozkładu LU
//---------------------------------------
bool lu(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;
}

// Program główny
//---------------

int main()
{
  double ** A, det;
  int n, i, j;
 
  cout << setprecision(4) << fixed;
 
  // odczytujemy stopień macierzy A
  cin >> n;

  // tworzymy macierz A o odpowiednich rozmiarach
  A = new double * [n];
  for(i = 0; i < n; i++)
    A[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];

  // obliczamy wyznacznik
  if(lu(n, A))
  {
    det = A[0][0];
    for(i = 1; i < n; i++)
      det *= A[i][i];
    cout << det << endl;
  }
  else
    cout << "DZIELNIK ZERO\n";

  // usuwamy macierz z pamięci

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

  return 0;
}

Basic

' Wyznacznik przez rozkład LU
' Data: 10.02.2011
' (C)2020 mgr Jerzy Wałaszek
'-----------------------------

Const eps As Double = 1e-12

' Funkcja realizuje algorytm rozkładu LU
'---------------------------------------
Function lu(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

  lu = 1

End Function

Dim As Integer i, j, n
Dim A As Double Ptr Ptr
Dim det As Double

Open Cons For Input As #1

' odczytujemy stopień macierzy A
Input #1, n

' tworzymy macierz A
' o odpowiednich rozmiarach
A = New Double Ptr[n]
For i = 0 To n-1
  A[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

' obliczamy wyznacznik
If lu(n, A) Then
  det = A[0][0]
  For i = 1 To n-1
    det *= A[i][i]
  Next
  Print det
Else
  Print "DZIELNIK ZERO"
End If

' usuwamy macierz z pamięci
For i = 0 To n-1
  Delete [] A[i]
Next
Delete [] A

End

Python (dodatek)

# Wyznacznik przez rozkład LU
# Data: 30.04.2024
# (C)2024 mgr Jerzy Wałaszek
# ----------------------------

eps = 1e-12


# Funkcja dokonuje rozkładu LU
# macierzy a
# -----------------------------
def lu(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


# Program główny
# ===============

# odczytujemy stopień macierzy a
n = int(input())
# odczytujemy macierz a
a = []
for i in range(n):
    s = input().split()
    arr = [float(s[j]) for j in range(n)]
    a.append(arr)
# obliczamy wyznacznik
if lu(a):
    det = a[0][0]
    for i in range(1, n):
        det *= a[i][i]
    print(det)
else:
    print("DZIELNIK ZERO")

Wynik:

Na początek: podrozdziału 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 – wyznacznik macierzy