Algorytmy i Struktury Danych - Rozkład LU

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 det lub informacja o błędzie, det ∈ R

Zmienne pomocnicze:

indeks, i ∈ N

Lista kroków:

K01:	Wyznacz rozkład LU macierzy A, wynik w A
K02:	Jeśli się powiodło, idź do kroku K05
K03:	Pisz "DZIELNIK ZERO"
K04:	Zakończ
K05:	det ← A [1, 1]	liczymy wyznacznik
K06:	Dla i = 2, 3, …, n : wykonuj det ← det × A [i, i ]	kolejne iloczyny wyrazów 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;

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

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

Wynik:

3
1 2 3
6 5 4
3 7 2
63.0000

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