Rozkład LU


Tematy pokrewne   Podrozdziały
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
  Rozwiązanie nr 1
Rozwiązanie nr 2
 

Problem

Daną macierz kwadratową An×n rozłożyć na dwie macierze: trójkątną dolną Ln×n o elementach głównej przekątnej równych 1 i trójkątną górną Un×n, takich że:

 

An×n = Ln×n × Un×n

 

 

 

Rozwiązanie nr 1

Rozważmy macierz kwadratową A o wymiarze n × n. Załóżmy, iż jest ona iloczynem dwóch macierzy trójkątnych L i U:

 

An × n =
  a1,1 a1,2 a1,3 ... a1,n-1 a1,n  
  a2,1 a2,2 a2,3 ... a2,n-1 a2,n  
  a3,1 a3,2 a3,3 ... a3,n-1 a3,n  
  ... ... ... ... ... ...  
  an-1,1 an-1,2 an-1,3 ... an-1,n-1 an-1,n  
  an,1 an,2 an,3 ... an,n-1 an,n  
= Ln × n × Un × n =
  1 0 0 ... 0 0  
  l2,1 1 0 ... 0 0  
  l3,1 l3,2 1 ... 0 0  
  ... ... ... ... ... ...  
  ln-1,1 ln-1,2 ln-1,3 ... 1 0  
  ln,1 ln,2 ln,3 ... ln,n-1 1  
×
  u1,1 u1,2 u1,3 ... u1,n-1 u1,n  
  0 u2,2 u2,3 ... u2,n-1 u2,n  
  0 0 u3,3 ... u3,n-1 u3,n  
  ... ... ... ... ... ...  
  0 0 0 ... un-1,n-1 un-1,n  
  0 0 0 ... 0 un,n  

 

Wykorzystując wzory na mnożenie macierzy otrzymujemy dla każdego elementu ai,j macierzy A następujące równanie:

 

 

Ponieważ macierze L i U są macierzami trójkątnymi, to część ich elementów jest równa zero. Otrzymujemy stąd dwa przypadki w zależności od indeksów i oraz j:

 

1)  i j, czyli element macierzy A ai,j leży na głównej przekątnej lub ponad nią. W takim przypadku ostatnim, niezerowym elementem w wierszu i-tym macierzy L jest element li,i. Skoro tak, to podany wzór możemy zredukować do postaci:

 

 

Dodatkowo z postaci macierzy L wiemy, iż element li,i, leżący na głównej przekątnej, jest zawsze równy 1. Stąd:

 

 

Wzór ten pozwala nam wyliczyć element ui,j o ile są znane poprzednie elementy macierzy L i U:

 

 

2) i > j, czyli element macierzy A ai,j leży pod główną przekątną. W tym przypadku ostatnim, niezerowym elementem w j-tej kolumnie macierzy U jest uj,j. Zatem wzór wyjściowy redukuje się do postaci:

 

 

Otrzymany wzór pozwala wyliczyć element li,j o ile znane są poprzednie elementy macierzy L i U:

 

 

W ten sposób doszliśmy do algorytmu Doolitle'a wyznaczania rozkładu macierzy A na iloczyn dwóch macierzy L i U (ang. LU matrix decomposition). Zasada działania tego algorytmu jest następująca:

Tworzymy dwie macierze L i U o stopniu macierzy A. Wypełniamy je zerami. W macierzy L główną przekątną wypełniamy elementami o wartości 1.

Dla j = 1,2,...,n wykonujemy kolejno dwie poniższe operacje:

        Dla i = 1,2,...,j wykorzystujemy wzór z punktu 1) do wyznaczenia elementów ui,j.

        Dla i = j+1, j+2,...,n wykorzystujemy wzór z punktu 2) do wyznaczenia elementów li,j.

 

Przykład:

W celu zobrazowania pracy algorytmu Doolitle'a wyznaczmy wzory na rozkład LU macierzy A3 × 3:

 

A3 × 3 =
  a1,1 a1,2 a1,3  
  a2,1 a2,2 a2,3  
  a3,1 a3,2 a3,3  
= L3 × 3 × U3 × 3 =
  1 0 0  
  l2,1 1 0  
  l3,1 l3,2 1  
×
  u1,1 u1,2 u1,3  
  0 u2,2 u2,3  
  0 0 u3,3  

 

j = 1

    i = 1, i = j, zatem wyznaczamy u1,1:
   

    i = 2, i > j, wyznaczamy l2,1:

   

    i = 3, i > j, wyznaczamy l3,1:

   

 

j = 2

    i = 1, i < j, wyznaczamy u1,2:

   

    i = 2, i = j, wyznaczamy u2,2:

   

    i = 3, i > j, wyznaczamy l3,2:

   

 

j = 3

    i = 1, i < j, wyznaczamy u1,3:

   

    i = 2, i < j, wyznaczamy u2,3:

   

    i = 3, i = j, wyznaczamy u3,3:

   

Dla przykładowych danych liczbowych otrzymujemy ze wzorów:

 

  1 4 7  
  2 5 8  
  3 6 9  
=
  1 0 0  
  2 1 0  
  3 2 1  
×
  1 4 7  
  0 -3 -6  
  0 0 0  

 

W powyższym przykładzie zwróć uwagę na następujące fakty:

Prowadzi to do wniosku, iż macierze L i U mogą zostać umieszczone w miejscu macierzy A, czyli operacja będzie wykonana in situ (łac. w miejscu) i zaoszczędzimy pamięć. W macierzy L istotna jest tylko część leżąca pod główną przekątną. Elementy przekątnej są zawsze równe 1, a elementy leżące nad przekątną są równe zero i nie musimy ich wcale zapamiętywać. Z drugiej strony w macierzy U istotna staje się tylko główna przekątna oraz elementy leżące nad nią – pozostałe elementy są równe zero i również nie musimy ich zapamiętywać. Zatem algorytm Doolitle'a może przekształcać macierz A na połączone macierze L i U:

 

  a1,1 a1,2 a1,3 ... a1,n-1 a1,n  
  a2,1 a2,2 a2,3 ... a2,n-1 a2,n  
  a3,1 a3,2 a3,3 ... a3,n-1 a3,n  
  ... ... ... ... ... ...  
  an-1,1 an-1,2 an-1,3 ... an-1,n-1 an-1,n  
  an,1 an,2 an,3 ... an,n-1 an,n  
  u1,1 u1,2 u1,3 ... u1,n-1 u1,n  
  l2,1 u2,2 u2,3 ... u2,n-1 u2,n  
  l3,1 l3,2 u3,3 ... u3,n-1 u3,n  
  ... ... ... ... ... ...  
  ln-1,1 ln-1,2 ln-1,3 ... un-1,n-1 un-1,n  
  ln,1 ln,2 ln,3 ... ln,n-1 un,n  

 

Algorytm Doolitle'a rozkładu LU w miejscu

Wejście
n    stopień macierzy, n N
A  – macierz stopnia n, elementy A R
Wyjście:

r = true, A – macierz zawiera połączone macierze L i U
r = false, błąd dzielenia przez zero.

Elementy pomocnicze:
i,j,k  –  indeksy i,j,k N
s  – suma iloczynów. s R
ε  – otoczenie zera, ε = 10-12, ε R
Lista kroków:
K01:  r ← false ; zakładamy porażkę
K02:  Dla j = 1,2,...,n wykonuj K03...K11  
K03:     Jeśli | A[j,j] | < ε, to zakończ ; rozkład niemożliwy, gdyż dzielnik równy zero
K04:     Dla i = 1,2,..., j wykonuj K05...K07 ; najpierw wyznaczamy w kolumnie wyrazy u
K05:         s ← 0  
K06:         Dla k = 1,2,...,i-1 wykonuj ss + A[i,k] × A[k,j] ; wyznaczamy sumę iloczynów wyrazów z L i U
K07:         A[i,j] ← A[i,j] - s ; obliczamy U[i,j]
K08:     Dla i = j+1, j+2,...,n wykonuj K09...K11  
K09:         s ← 0  
K10:         Dla k = 1,2,...,j-1 wykonuj ss + A[i,k] × A[k,j] ; wyznaczamy sumę iloczynów wyrazów z L i U
K11:
        A[i,j] ←   A[i,j] - s  
A[j,j]
; obliczamy L[i,j]
K12: r ← true ; zaznaczamy sukces
K13: Zakończ  

 

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 odczytuje ze standardowego wejścia macierz A (najpierw liczba n określająca stopień macierzy, a następnie n2 liczb dla kolejnych elementów),  wyznacza jej rozkład LU za pomocą algorytmu Doolitle'a i wyświetla wynikową macierz A. Przykładowe dane dla programu:

 

3
1 4 7
2 5 8
3 6 9

 

Lazarus
// Rozkład LU algorytmem Doolitle'a
// Data: 21.03.2010
// (C)2012 mgr Jerzy Wałaszek
//-----------------------------

program ludoolitle;

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

const

  eps = 1e-12;

// Funkcja realizuje algorytm Doolitle'a rozkładu LU
//--------------------------------------------------
function Doolitle(n : integer; A : MReal) : boolean;

var
  i,j,k : integer;
  s     : double;

begin

  for j := 0 to n - 1 do
  begin
    if abs(A[j][j]) < eps then Exit(false);
    for i := 0 to j do
    begin
      s := 0;
      for k := 0 to i - 1 do
        s := s + A[i][k] * A[k][j];
      A[i][j] := A[i][j] - s;
    end;
    for i := j + 1 to n - 1 do
    begin
      s := 0;
      for k := 0 to j - 1 do
        s := s + A[i][k] * A[k][j];
      A[i][j] := (A[i][j] - s) / A[j][j];
    end;
  end;
  Doolitle := true;
end;

// Program główny

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

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]);

  // wyznaczamy rozkład LU macierzy A

  if Doolitle(n,A) then
    for i := 0 to n - 1 do
    begin
      for j := 0 to n - 1 do write(A[i][j]:8:3,' ');
      writeln;
    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.
Code::Blocks
// Rozkład LU algorytmem Doolitle'a
// Data: 21.03.2010
// (C)2012 mgr Jerzy Wałaszek
//-----------------------------

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

using namespace std;

const double eps = 1e-12;

// Funkcja realizuje algorytm Doolitle'a rozkładu LU
//--------------------------------------------------
bool Doolitle(int n, double ** A)
{
  int i,j,k;
  double s;
  
  for(j = 0; j < n; j++)
  {
    if(fabs(A[j][j]) < eps) return false;
    for(i = 0; i <= j; i++)
    {
      s = 0;
      for(k = 0; k < i; k++) s +=  A[i][k] * A[k][j];
      A[i][j] -= s;
    }
    for(i = j + 1; i < n; i++)
    {
      s = 0;
      for(k = 0; k < j; k++) s += A[i][k] * A[k][j];
      A[i][j] = (A[i][j] - s) / A[j][j];
    }
  }
  return true;
}

// Program główny

int main()
{
  double **A;
  int n,i,j;
  
  cout << setprecision(3) << 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];

  // wyznaczamy rozkład LU macierzy A

  if(Doolitle(n,A))
    for(i = 0; i < n; i++)
    {
      for(j = 0; j < n; j++)
        cout << setw(8) << A[i][j] << " ";
      cout << endl;
    }
  else
    cout << "DZIELNIK ZERO\n";

  // usuwamy macierz z pamięci

  for(i = 0; i < n; i++) delete [] A[i];
  delete [] A;
  return 0;
}
Free Basic
' Rozkład LU algorytmem Doolitle'a
' Data: 21.03.2010
' (C)2012 mgr Jerzy Wałaszek
'-----------------------------

Const eps As Double = 1e-12

' Funkcja realizuje algorytm Doolitle'a rozkładu LU
'--------------------------------------------------
Function doolitle(n As Integer, A As Double Ptr Ptr) As Integer

  Dim As Integer i,j,k
  Dim As Double s
  
  For j = 0 To n - 1
    If Abs(A[j][j]) < eps Then Return 0
    For i = 0 To j
      s = 0
      For k = 0 To i - 1
        s += A[i][k] * A[k][j]
      Next
      A[i][j] -= s
    Next
    For i = j + 1 To n - 1
      s = 0
      For k = 0 To j - 1
        s += A[i][k] * A[k][j]
      Next
      A[i][j] = (A[i][j] - s) / A[j][j]
    Next
  Next
  doolitle = 1
End Function

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

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

' wyznaczamy rozkład LU macierzy A

If doolitle(n,A) = 1 Then
  For i = 0 To n - 1
    For j = 0 To n - 1
      Print Using "####.### ";A[i][j];
    Next
    Print
  Next
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 4 7
2 5 8
3 6 9
   1.000    4.000    7.000
   2.000   -3.000   -6.000
   3.000    2.000    0.000

 

Rozwiązanie nr 2

W literaturze informatycznej można znaleźć prostsze rozwiązanie wyznaczania rozkładu LU macierzy. Polega ono na przechodzeniu przez kolejne elementy głównej przekątnej macierzy A od pierwszego do przedostatniego i dla każdego elementu A[k,k] wykonywaniu kolejno dwóch operacji:

  1. Normalizacji kolumny – elementy w kolumnie pod elementem A[k,k] są dzielone przez ten element.
  2. Modyfikacji podmacierzy, która zawiera wszystkie elementy macierzy A leżące poniżej i na prawo elementu przekątnej z pominięciem wiersza i kolumny zawierającego element przekątnej. Modyfikacja polega na odjęciu od każdego elementu A[i.j] tej podmacierzy iloczynu elementów A[i,k] × A[k,j].

Dzięki takiemu podejściu algorytm rozkładu LU znacznie się upraszcza.

 

Algorytm rozkładu LU w miejscu

Wejście
n    stopień macierzy, n N
A  – macierz stopnia n, elementy A R
Wyjście:

r = true, A – macierz zawiera połączone macierze L i U
r = false, błąd dzielenia przez zero.

Elementy pomocnicze:
i,j,k  –  indeksy i,j,k N
ε  – otoczenie zera, ε = 10-12, ε R
Lista kroków:
K01:  r ← false ; zakładamy porażkę
K02:  Dla k = 1,2,...,n - 1 wykonuj K03...K06  
K03:     Jeśli | A[k,k] | < ε, to zakończ ; element zerowy nie pozwoli wykonać dzielenia
K04:     Dla i = k+1,k+2,...,n wykonuj:
     A[i,j] ←   A[i,k]
A[k,k]
; normalizujemy kolumnę
K05:     Dla i = k+1,k+2,...,n wykonuj K06  
K06:         Dla j = k+1,k+2,...,n wykonuj:
            A[i,j] ← A[i,j] - A[i,k] × A[k,j]
; modyfikujemy podmacierz
K07:  r ← true ; sukces
K08: Zakończ  

 

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 odczytuje ze standardowego wejścia macierz A (najpierw liczba n określająca stopień macierzy, a następnie n2 liczb dla kolejnych elementów),  wyznacza jej rozkład LU  i wyświetla wynikową macierz A.

 

Lazarus
// Rozkład LU
// Data: 30.04.2010
// (C)2012 mgr Jerzy Wałaszek
//-----------------------------

program ludecomp;

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;

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]);

  // wyznaczamy rozkład LU macierzy A

  if lu(n,A) then
    for i := 0 to n - 1 do
    begin
      for j := 0 to n - 1 do write(A[i][j]:8:3,' ');
      writeln;
    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.
Code::Blocks
// Rozkład LU
// Data: 30.04.2010
// (C)2012 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;
  int n,i,j;
  
  cout << setprecision(3) << 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];

  // wyznaczamy rozkład LU macierzy A

  if(lu(n,A))
    for(i = 0; i < n; i++)
    {
      for(j = 0; j < n; j++)
        cout << setw(8) << A[i][j] << " ";
      cout << endl;
    }
  else
    cout << "DZIELNIK ZERO\n";

  // usuwamy macierz z pamięci

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

  return 0;
}
Free Basic
' Rozkład LU
' Data: 30.04.2010
' (C)2012 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 A As Double Ptr Ptr
Dim As Integer i,j,n

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

' wyznaczamy rozkład LU macierzy A

If lu(n,A) = 1 Then
  For i = 0 To n - 1
    For j = 0 To n - 1
      Print Using "####.### ";A[i][j];
    Next
    Print
  Next
Else
    Print "DZIELNIK ZERO"
End If

' Usuwamy macierz A

For i = 0 To n - 1
  Delete [] A[i]
Next
Delete [] A

End

 



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.