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


Tematy pokrewne
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
 

Problem

Znaleźć rozwiązanie układu równań liniowych postaci:

 

a1,1x1 + a1,2x2 + ... + a1,nxn = b1
a2,1x1 + a2,2x2 + ... + a2,nxn = b2
...
an,1x1 + an,2x2 + ... + an,nxn = bn

 

wykorzystując rozkład LU.

 

 

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

 

 

lub w skrócie:

 

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  
  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  

 

U =
  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  

 

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):


 

podstawiając wstecz (ang. back substitution):

 

 

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 K04...K06  
K04:     s ← 0 ; zerujemy sumę
K05:     Dla j = 1,2,...,i-1: wykonuj
        ss + A[i,j] × X[j]
; obliczamy sumę iloczynów l(i,j) x y(j)
K06:     X[i] ← B[i] - s ; obliczamy y(i)
K07: Jeśli |A[n,n]| < ε, to zakończ ; błąd dzielenia przez zero
K08:
X[n] ←   X[n]
A[n,n]
; teraz obliczamy wektor X
K09: Dla i = n-1,n-2,...,1: wykonuj K10...K13  
K10:     s ← 0 ; zerujemy sumę
K11:     Dla j = i+1,i+2,...,n: wykonuj
        ss + A[i,j] × X[j]
; obliczamy sumę iloczynów u(i,j) x x(j)
K12:     Jeśli |A[i,i]| < ε, to zakończ ; błąd dzielenia przez zero
K13:
    X[i] ←   X[i] - s
A[i,i]
; obliczamy x(i)
K14:  r ← true ; zaznaczamy sukces
K15: 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 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:

 

4x1
3x1
2x1
2x1
 - 2x2
+   x2
 + 4x2
 - 2x2
 + 4x3
+ 4x3
 + 2x3
+ 4x3
 - 2x4
+ 2x4
 + x4
+ 2x4
 =  8
=  7
=10
=  2

 

Dane dla programu:

 

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

 

Lazarus
// Rozkład LU - układ równań liniowych
// Data: 18.10.2010
// (C)2012 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.
Code::Blocks
// Rozkład LU - układ równań liniowych
// Data: 18.10.2010
// (C)2012 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;
}
Free Basic
' Rozkład LU - układ równań liniowych
' Data: 18.10.2010
' (C)2012 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
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
 


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.