Algorytmy i Struktury Danych

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

A_n×n =

a_1,1	a_1,2	a_1,3	…	a_1,n-1	a_1,n
a_2,1	a_2,2	a_2,3	…	a_2,n-1	a_2,n
a_3,1	a_3,2	a_3,3	…	a_3,n-1	a_3,n
…	…	…	…	…	…
a_n-1,1	a_n-1,2	a_n-1,3	…	a_n-1,n-1	a_n-1,n
a_n,1	a_n,2	a_n,3	…	a_n,n-1	a_n,n

= L_n×n×U_n×n =

1	0	0	…	0	0
l_2,1	1	0	…	0	0
l_3,1	l_3,2	1	…	0	0
…	…	…	…	…	…
l_n-1,1	l_n-1,2	l_n-1,3	…	1	0
l_n,1	l_n,2	l_n,3	…	l_n,n-1	1

×

u_1,1	u_1,2	u_1,3	…	u_1,n-1	u_1,n
0	u_2,2	u_2,3	…	u_2,n-1	u_2,n
0	0	u_3,3	…	u_3,n-1	u_3,n
…	…	…	…	…	…
0	0	0	…	u_n-1,n-1	u_n-1,n
0	0	0	…	0	u_n,n

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

A_n×n= L_n×n×U_n×n

a_i,j =

^k=1

l_i,k×u_k,j

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 a_i,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 l_i,i. Skoro tak, to podany wzór możemy zredukować do postaci:

a_i,j =

^k=1

l_i,k×u_k,j

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

a_i,j = (

_i-1

^k=1

l_i,k×u_k,j)+u_i,j

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

u_i,j = a_i,j-

_i-1

^k=1

l_i,k×u_k,j

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

a_i,j =

^k=1

l_i,k×u_k,j

a_i,j = (

_j_-1

^k=1

l_i,k×u_k,j)+l_i,j×u_{j, j}

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

l_i,j =

×(a_i,j-

_i-1

Σ

^k=1

l_i,k×u_k,j)

u_j,j

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 u_{i, j}.

Dla i = j+1, j+2, …, n: wykorzystujemy wzór z punktu 2) do wyznaczenia elementów l_{i, j}.

Przykład:

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

A_3×3 =

a_1,1a_1,2a_1,3

a_2,1a_2,2a_2,3

a_3,1a_3,2a_3,3

= L_3×3×U_3×3=

1	0	0
l_2,1	1	0
l_3,1	l_3,2	1

×

u_1,1	u_1,2	u_1,3
0	u_2,2	u_2,3
0	0	u_3,3

j = 1

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

u_1,1 = a_1,1-

₀

^k=1

l_1,k×u_k,1 = a_1,1

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

l_2,1 = 1/u_1,1×(a_2,1-

₀

^k=1

l_2,k×u_k,1) = 1/u_1,1×a_2,1

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

l_3,1 = 1/u_1,1×(a_3,1-	₀	l_3,k×u_k,1) = 1/u_1,1×a_3,1
	Σ
	^k=1

j = 2

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

u_1,2 = a_1,2-

₀

^k=1

l_1,k×u_k,2 = a_1,2

i = 2, i = j, wyznaczamy u_2,2:

u_2,2 = a_2,2-

₁

^k=1

l_2,k×u_k,2 = a_2,2-l_2,1×u_1,2

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

l_3,2 = 1/u_2,2×(a_3,2-

₁

^k=1

l_3,k×u_k,2) = 1/u_2,2×(a_3,2-l_3,1×u_1,2)

j = 3

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

u_1,3 = a_1,3-

₀

^k=1

l_1,k×u_k,3 = a_1,3

i = 2, i < j, wyznaczamy u_{2, 3}:

u_2,3 = a_2,3-

₁

^k=1

l_2,k×u_k,3 = a_2,3-l_2,1×u_1,3

i = 3, i = j, wyznaczamy u_3,3:

u_3,3 = a_3,3-

₂

^k=1

l_3,k×u_k,3 = a_3,3-l_3,1×u_1,3-l_3,2×u_2,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:

Potrzebne we wzorach elementy macierzy L i U są zawsze wyliczone wcześniej.
We wszystkich wzorach elementy macierzy A są używane tylko jeden raz.

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:

a_1,1	a_1,2	a_1,3	…	a_1,n-1	a_1,n
a_2,1	a_2,2	a_2,3	…	a_2,n-1	a_2,n
a_3,1	a_3,2	a_3,3	…	a_3,n-1	a_3,n
…	…	…	…	…	…
a_n-1,1	a_n-1,2	a_n-1,3	…	a_n-1,n-1	a_n-1,n
a_n,1	a_n,2	a_n,3	…	a_n,n-1	a_n,n

→

→

u_1,1	u_1,2	u_1,3	…	u_1,n-1	u_1,n
l_2,1	u_2,2	u_2,3	…	u_2,n-1	u_2,n
l_3,1	l_3,2	u_3,3	…	u_3,n-1	u_3,n
…	…	…	…	…	…
l_n-1,1	l_n-1,2	l_n-1,3	…	u_n-1,n-1	u_n-1,n
l_n,1	l_n,2	l_n,3	…	l_n,n-1	u_n,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, macierz A 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 kroki K03…K11
K03:   Jeśli |A[j,j]| < ε, ; rozkład niemożliwy, 
       to zakończ ; gdyż dzielnik równy zero
K04:   Dla i = 1, 2, …, j: ; najpierw wyznaczamy w kolumnie wyrazy u
       wykonuj kroki K05…K07
K05:     s ← 0
K06:     Dla k = 1, 2, …, i-1: ; wyznaczamy sumę iloczynów wyrazów z L i U
         wykonuj s ← s+A[i,k]×A[k,j]
K07:     A[i,j] ← A[i,j]-s ; obliczamy U[i,j]
K08:   Dla i = j+1, j+2, …, n:
       wykonuj kroki K09…K11
K09:     s ← 0
K10:     Dla k = 1, 2, …, j-1: ; wyznaczamy sumę iloczynów wyrazów z L i U
         wykonuj s ← s+A[i,k]×A[k,j]
K11:     A[i,j] ← (A[i,j]-s):A[j,j] ; obliczamy L[i,j]
K12: r ← true ; zaznaczamy sukces
K13: 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 wyświetla wynikową macierz A.

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

3
1 4 7
2 5 8
3 6 9

Pascal

// Rozkład LU algorytmem Doolitle'a
// Data: 21.03.2010
// (C)2020 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.

C++

// Rozkład LU algorytmem Doolitle'a
// Data: 21.03.2010
// (C)2020 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;
}

Basic

' Rozkład LU algorytmem Doolitle'a
' Data: 21.03.2010
' (C)2020 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

Python (dodatek)

# Rozkład LU algorytmem Doolitle'a
# Data: 23.04.2024
# (C)2024 mgr Jerzy Wałaszek
#-----------------------------

EPS = 1e-12


# Funkcja realizuje algorytm
# Doolitle'a rozkładu LU
#---------------------------
def doolitle(a):
    n = len(a)
    for j in range(n):
        if abs(a[j][j]) < EPS:
            return False
        for i in range(j+1):
            s = 0
            for k in range(i):
                s += a[i][k]*a[k][j]
            a[i][j] -= s
        for i in range(j+1, n):
            s = 0
            for k in range(j):
                s += a[i][k]*a[k][j]
            a[i][j] = (a[i][j]-s)/a[j][j]
    return True


# Program główny

# odczytujemy stopień macierzy A
n = int(input())
# odczytujemy macierz A
a = []
for i in range(n):
    arr = input().split()
    arr = [float(arr[j]) for j in range(n)]
    a.append(arr)
# wyznaczamy rozkład LU macierzy A
if doolitle(a):
    for i in range(n):
        for j in range(n):
            print("%8.3f " % a[i][j], end="")
        print()
else:
    print("DZIELNIK ZERO")
# usuwamy macierz z pamięci
a = None

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

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

Problem

Rozwiązanie nr 1

j = 1

j = 2

j = 3

Algorytm Doolitle'a rozkładu LU w miejscu

Wejście:

Wyjście:

Elementy pomocnicze:

Lista kroków:

Przykładowe programy

Rozwiązanie nr 2

Algorytm rozkładu LU w miejscu

Wejście:

Wyjście:

Elementy pomocnicze:

Lista kroków:

Przykładowe programy