Serwis Edukacyjny
w I-LO w Tarnowie
obrazek

Materiały dla uczniów liceum

  Wyjście       Spis treści       Wstecz       Dalej  

Autor artykułu: mgr Jerzy Wałaszek

©2020 mgr Jerzy Wałaszek
I LO w Tarnowie

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

SPIS TREŚCI
Tematy pomocnicze

Problem

Należy znaleźć rozwiązanie układu równań liniowych postaci:

wykorzystując rozkład LU

Rozwiązanie

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

obrazek

podstawiając wstecz ( ang. back substitution ):

obrazek

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, nN
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.

Zmienne pomocnicze:

i, j,  –  indeksy, i, jN
ε  – otoczenie zera, ε  = 10 -12, ε  ∈ R
s  – przechowuje sumę iloczynów, sR

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:
    wykonuj s  ← s  + 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: 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 :
    wykonuj s  ← s  + A [ i, j  ] × X [ j  ]
obliczamy sumę iloczynów u ( i, j ) · x ( j )
K12:     Jeśli | A [ i, i  ]  | < ε,
    to zakończ
błąd dzielenia przez zero
K13:     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 1
3x 1
2x 1
2x 1
 - 2x 2
+   x 2
 + 4x 2
 - 2x 2
 + 4x 3
+ 4x 3
 + 2x 3
+ 4x 3
 - 2x 4
+ 2x 4
 + x 4
+ 2x 4
=  8
=  7
=10
=  2
Dane przykładowe ( przekopiuj do schowka i wklej do okna konsoli ):

4
4 -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
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
Na początek:  podrozdziału   strony 

Zespół Przedmiotowy
Chemii-Fizyki-Informatyki

w I Liceum Ogólnokształcącym
im. Kazimierza Brodzińskiego
w Tarnowie
ul. Piłsudskiego 4
©2020 mgr Jerzy Wałaszek

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.