Algorytmy i Struktury Danych - Rozkład LU

Układ równań liniowych możemy zapisać w postaci macierzowej jako:

A×X = B

gdzie:

A – macierz kwadratowa o wymiarze n×n współczynników stojących przy niewiadomych x.
X – wektor kolumnowy n niewiadomych.
B – wektor kolumnowy n wyrazów wolnych.

Załóżmy, że macierz współczynników A możemy przedstawić jako iloczyn dwóch macierzy L i U:

A = L×U

gdzie:

L =

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 =

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

Wtedy:

A×X = B
(L×U)×X = B
L×(U×X) = B

Najpierw znajdziemy wektor Y, taki że:

L×Y = B

A następnie mając Y obliczymy wektor X, który jest rozwiązaniem układu równań:

U×X = Y

Korzyścią tego rozwiązania jest to, iż macierze L i U ze względu na swoją postać pozwalają na bardzo proste wyliczenie Y i X.

podstawiając w przód (ang. forward substitution):

y₁= b₁

y_i = b_i -

^i-1

^j=1

l_i,j×y_j, dla i = 2, 3, …, n

podstawiając wstecz (ang. back substitution):

x_n=	y_n
	u_n,n

y_i =

×(y_i-

ⁿ

^j=i+1

u_i,j×x_j), dla i = n-1, n-2, …, 1

u_i,i

Wektor pomocniczy Y może być w pamięci komputera przechowywany w wektorze X, co zaoszczędzi miejsca.

To podejście szczególnie jest przydatne w przypadku konieczności rozwiązywania wielu układów równań, które różnią się jedynie wektorem wyrazów wolnych B. W takim przypadku dokonujemy jednokrotnie rozkładu LU macierzy współczynników A, a następnie wymieniając jedynie wektor B rozwiązujemy kolejne układy przy pomocy wyliczonych wcześniej macierzy L i U.

Algorytm rozwiązywania układu równań liniowych po rozkładzie LU

Wejście:

n : stopień macierzy; n ∈ N.
A : macierz kwadratowa stopnia n z połączonymi macierzami L i U otrzymana w algorytmie rozkładu LU; Elementy A ∈ R.
B : wektor stopnia n, zawierający wyrazy wolne układu równań; Elementy B ∈ R.

Wyjście:

r = true, wektor X zawiera rozwiązanie układu równań.
r = false, błąd dzielenia przez zero.

Elementy pomocnicze:

i, j, : indeksy; i, j ∈ N.
ε : otoczenie zera; ε = 10^-12, ε ∈ R.
s : przechowuje sumę iloczynów; s ∈ R.

Lista kroków:

K01: r ← false ; zakładamy porażkę
K02: X[1] ← B[1] ; najpierw wyliczamy wyrazy wektora Y
K03: Dla i = 2, 3, …, n:
     wykonuj kroki K04…K06
K04:   s ← 0 ; zerujemy sumę
K05:   Dla j = 1, 2, …, i-1: ; obliczamy sumę
       wykonuj s ← s+A[i, j]×X[j] ; iloczynów l[i, j]×y[j]
K06:   X[i] ← B[i]-s ; obliczamy y[i]
K07: Jeśli |A[n, n]| < ε, ; błąd dzielenia przez zero
     to zakończ
K08: X[i] ← X[n]:A[n, n] ; teraz obliczamy wektor X
K09: Dla i = n-1, n-2, …, 1:
     wykonuj kroki K10…K13
K10:   s ← 0 ; zerujemy sumę
K11:   Dla j = i+1, i+2, …, n: ; obliczamy sumę
       wykonuj s ← s+A[i, j]×X[j] ; iloczynów u[i, j]×x[j]
K12:   Jeśli |A[i, i]| < ε,  ; błąd dzielenia przez zero
       to zakończ
K13:   X[i] ← (X[i]-s):A[i, i] ; obliczamy x[i]
K14: r ← true ; zaznaczamy sukces
K15: 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 wczytuje dane definiujące macierz współczynników AB i wylicza niewiadome X. Jeśli obliczenie nie może zostać wykonane, program wypisuje komunikat "DZIELNIK ZERO".

Dane dla macierzy AB są zdefiniowane następująco:

Pierwsza liczba określa liczbę niewiadomych n.

Następne n×(n+1) liczb określa zawartość macierzy AB wierszami – pierwsze n liczb każdego wiersza odnosi się do współczynników a, a (n+1)-sza liczba odnosi się do wyrazu wolnego b. Poniżej mamy przykładowe dane wejściowe dla programu:

Równanie:

4x₁–2x₂+4x₃–2x₄=  83x₁+1x₂+4x₃+2x₄=  72x₁+4x₂+2x₃+1x₄= 102x₁–2x₂+4x₃+2x₄=  2

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

44 -2 4 -2 8
3 1 4 2 7
2 4 2 1 10
2 -2 4 2 2

Pascal

// Rozkład LU-układ równań liniowych
// Data: 18.10.2010
// (C)2020 mgr Jerzy Wałaszek
//------------------------------------

program lueq;

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

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(n : integer;
                 A : MReal;
                 B : NReal;
                 X : NReal)
                 : boolean;
var
  i, j : integer;
  s   : double;
begin

  X[0] := B[0];

  for i := 1 to n-1 do
  begin
    s := 0;
    for j := 0 to i-1 do
      s := s+A[i][j]*X[j];
    X[i] := B[i]-s;
  end;

  if abs(A[n-1][n-1]) < eps then
    Exit(false);

  X[n-1] := X[n-1]/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];
    if abs(A[i][i]) < eps then
      Exit(false);
    X[i] := (X[i]-s)/A[i][i];
  end;

  lusolve := true;

end;

// Program główny

var
  A     : MReal;
  B, X   : NReal;
  n, i, j : integer;

begin

  // odczytujemy liczbę niewiadomych
  read(n);

  // tworzymy macierze A, B i X
  SetLength(A, n);
  SetLength(B, n);
  SetLength(X, n);

  for i := 0 to n-1 do
    SetLength(A[i], n);

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

  // rozwiązujemy układ i wyświetlamy wyniki
  if ludist(n, A) and lusolve(n, A, B, X) then
  begin
    for i := 0 to n-1 do
      writeln('x', i+1, ' = ', X[i]:9:4);
  end
  else
    writeln('DZIELNIK ZERO');

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

end.

C++

// Rozkład LU-układ równań liniowych
// Data: 18.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 B
//--------------------------
bool lusolve(int n, 
             double ** A, 
             double * B, 
             double * X)
{
  int i, j;
  double s;

  X[0] = B[0];

  for(i = 1; i < n; i++)
  {
    s = 0;
    for(j = 0; j < i; j++)
      s += A[i][j]*X[j];
    X[i] = B[i]-s;
  }

  if(fabs(A[n-1][n-1]) < eps)
    return false;

  X[n-1] /= 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];
    if(fabs(A[i][i]) < eps)
      return false;
    X[i] = (X[i]-s)/A[i][i];
  }

  return true;
}

// Program główny

int main()
{
  double **A, *B, *X;
  int n, i, j;

  cout << setprecision(4) << fixed;

  // odczytujemy liczbę niewiadomych
  cin >> n;

  // tworzymy macierze A, B i X
  A = new double * [n];
  B = new double[n];
  X = new double[n];
  for(i = 0; i < n; i++)
    A[i] = new double[n];

  // odczytujemy dane dla macierzy A i B
  for(i = 0; i < n; i++)
  {
    for(j = 0; j < n; j++)
      cin >> A[i][j];
    cin >> B[i];
  }

  // rozwiązujemy układ i wyświetlamy wyniki
  if(ludist(n, A) && lusolve(n, A, B, X))
  {
    for(i = 0; i < n; i++)
      cout << "x" << i+1 << " = "
           << setw(9) << X[i] << endl;
  }
  else
    cout << "DZIELNIK ZERO\n";

  // usuwamy macierze z pamięci
  for(i = 0; i < n; i++)
    delete [] A[i];
  delete [] A;
  delete [] B;
  delete [] X;

  return 0;
}

Basic

' Rozkład LU-układ równań liniowych
' Data: 18.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 B
'--------------------------
Function lusolve (n As Integer, _
                  A As Double Ptr Ptr, _
                  B As Double Ptr, _
                  X As Double Ptr) _
                  As Integer

  Dim As Integer i, j
  Dim As Double s

  X[0] = B[0]

  For i = 1 To n-1
    s = 0

    For j = 0 To i-1
      s += A[i][j]*X[j]
    Next

    X[i] = B[i]-s
  Next

  If Abs(A[n-1][n-1]) < eps Then _
    Return 0

  X[n-1] /= 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]
    Next

    If Abs(A[i][i]) < eps Then _
      Return 0

    X[i] = (X[i]-s)/A[i][i]
  Next

  lusolve = 1

End Function

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

Open Cons For Input As #1

' odczytujemy liczbę niewiadomych
Input #1, n

' tworzymy macierze A, B i X
A = New Double Ptr[n]
B = New Double[n]
X = New Double[n]
 
For i = 0 To n-1
  A[i] = New Double[n]
Next

' odczytujemy dane dla macierzy A i B
For i = 0 To n-1
  For j = 0 To n-1
    Input #1, A[i][j]
  Next
  Input #1, B[i]
Next

Close #1

' rozwiązujemy układ i wyświetlamy wyniki
If ludist(n, A) And lusolve(n, A, B, X) Then
  For i = 0 To n-1
    Print Using "x# = ####.####";i+1;X[i]
  Next
Else
  Print "DZIELNIK ZERO"
End If

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

End

Python (dodatek)

# Rozkład LU-układ równań liniowych
# Data: 25.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 b
#--------------------------
def lusolve(a, b, x):
    n = len(a)
    x[0] = b[0]
    for i in range(1, n):
        s = 0
        for j in range(i):
            s += a[i][j]*x[j]
        x[i] = b[i]-s
    if abs(a[n-1][n-1]) < EPS:
        return False
    x[n-1] /= a[n-1][n-1]
    for i in reversed(range(n-1)):
        s = 0
        for j in range(i+1, n):
            s += a[i][j]*x[j]
        if abs(a[i][i]) < EPS:
            return False
        x[i] = (x[i]-s)/a[i][i]
    return True


# Program główny
#---------------

# odczytujemy liczbę niewiadomych
n = int(input())
# tworzymy macierze a, b i x
a = []
b = []
x = [0]*n
for i in range(n):
    s = input().split()
    arr = [float(s[j]) for j in range(n+1)]
    b.append(arr.pop())
    a.append(arr)
# rozwiązujemy układ i wyświetlamy wyniki
if ludist(a) and lusolve(a, b, x):
    for i in range(n):
        print("x%d = %9.4f" % (i+1, x[i]))
else:
    print("DZIELNIK ZERO")

Wynik:

4
4 -2 4 -2 8
3 1 4 2 7
2 4 2 1 10
2 -2 4 2 2
x1 =   -1.0000
x2 =    2.0000
x3 =    3.0000
x4 =   -2.0000

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 – rozwiązywanie układu równań liniowych