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

©2024 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:

a1,1·x1+a1,2·x2+…+a1,n·xn = b1
a2,1·x1+a2,2·x2+…+a2,n·xn = b2
            
an,1·x1+an,2·x2+…+an,n·xn = bn

wykorzystując rozkład LU

Rozwiązanie

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 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×UX = 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 LU ze względu na swoją postać pozwalają na bardzo proste wyliczenie YX.

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

y1 = b1
yi = bi -

i-1

Σ

j=1

li,j×yj, dla i=2,3,…,n

podstawiając wstecz (ang. back substitution):

xn = 
 yn
un,n
yi = 
 1
×(yi -

n

Σ 

j=i+1

ui,j×xj), dla i=n-1,n-2,…,1
ui,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 LU.

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 LU 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ń.
= 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ę iloczynów l[i,j]×y[j]
       wykonuj s ← s+A[i,jX[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ę iloczynów u[i,j]×x[j]
       wykonuj s ← s+A[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:

4x1–2x2+4x3–2x4 =  8
3x1+1x2+4x3+2x4 =  7
2x1+4x2+2x3+1x4 = 10
2x1–2x2+4x3+2x4 =  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
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):
    a = input().split()
    a = [float(a[j]) for j in range(n+1)]
    B.append(a.pop())
    A.append(a)
    del a


# 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")

# usuwamy macierze z pamięci

del A,B,X
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
©2024 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.