Algorytmy i Struktury Danych - Rozkład LU

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

\(A\times X = B\)

gdzie:

\(A\) – macierz kwadratowa o wymiarze \({n\times 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\times U\)

gdzie:

\[
L =
\begin{bmatrix}
1 & 0 & 0 & \dots & 0 & 0 \\
l_{2,1} & 1 & 0 & \dots & 0 & 0 \\
l_{3,1} & l_{3,2} & 1 & \dots & 0 & 0 \\
\dots & \dots & \dots & \dots & \dots & \dots \\
l_{n-1,1} & l_{n-1,2} & l_{n-1,3} & \dots & 1 & 0 \\
l_{n,1} & l_{n,2} & l_{n,3} & \dots & l_{n,n-1} & 1
\end{bmatrix}
\\ \\
U =
\begin{bmatrix}
u_{1,1} & u_{1,2} & u_{1,3} & \dots & u_{1,n-1} & u_{1,n} \\
0 & u_{2,2} & u_{2,3} & \dots & u_{2,n-1} & u_{2,n} \\
0 & 0 & u_{3,3} & \dots & u_{3,n-1} & u_{3,n} \\
\dots & \dots & \dots & \dots & \dots & \dots \\
0 & 0 & 0 & \dots & u_{n-1,n-1} & u_{n-1,n} \\
0 & 0 & 0 & \dots & 0 & u_{n,n}
\end{bmatrix}
\]

Wtedy:

\[
\begin{aligned}
&A \times X = B \\
&(L \times U) \times X = B \\
&L \times (U \times X) = B
\end{aligned}
\]

Najpierw znajdziemy wektor \(Y\), taki że:

\[
L\times Y = B
\]

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

\[
U\times 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):

\[
\begin{aligned}
&y_{1} = b_{1} \\
&y_{i} = b_{i} - \sum\limits_{j=1}^{i-1} l_{i,j} \cdot y_{j}, \text{ dla } i = 2, 3, \dots, n
\end{aligned}
\]

podstawiając wstecz (ang. back substitution):

\[
\begin{aligned}
&x_{n} = \frac{y_{n}}{u_{n,n}} \\
&x_{i} = \frac{1}{u_{i,i}} \cdot \left( y_{i} - \sum\limits_{j=i+1}^{n} u_{i,j} \cdot x_{j} \right), \text{ dla } i = n-1, n-2, \dots \rightarrow 1
\end{aligned}
\]

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\in\mathbb{N}}\).
\(A\) : macierz kwadratowa stopnia \(n\) z połączonymi macierzami \(L\) i \(U\) otrzymana w algorytmie rozkładu LU; \({elementy \; A\in\mathbb{R}}\).
\(B\) : wektor stopnia \(n\), zawierający wyrazy wolne układu równań; \({elementy \; B\in\mathbb{R}}\).

Wyjście:

\(r = \text{true}\), wektor \(X\) zawiera rozwiązanie układu równań.
\(r = \text{false}\), błąd dzielenia przez zero.

Elementy pomocnicze:

\(i, j\) : indeksy; \({i,j\in\mathbb{N}}\).
\(\varepsilon\) : otoczenie zera; \({\varepsilon = 10^{-12}, \; \varepsilon\in\mathbb{R}}\).
\(s\) : przechowuje sumę iloczynów; \({s\in\mathbb{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 ilość niewiadomych \(n\).

Następne \({n\cdot(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){-}\text{sza}}\) liczba odnosi się do wyrazu wolnego \(b\). Poniżej mamy przykładowe dane wejściowe dla programu:

Równanie:

\[
\begin{cases}
4x_{1} - 2x_{2} + 4x_{3} - 2x_{4} = \phantom{1}8 \\
3x_{1} + 1x_{2} + 4x_{3} + 2x_{4} = \phantom{1}7 \\
2x_{1} + 4x_{2} + 2x_{3} + 1x_{4} = 10 \\
2x_{1} - 2x_{2} + 4x_{3} + 2x_{4} = \phantom{1}2
\end{cases}
\]

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

Rozkład LU – rozwiązywanie układu równań liniowych