KATEDRA URZ
DZEC
I SYSTEMÓW AUTOMATYKI
LABORATORIUM
METOD NUMERYCZYCH
W TECHNICE
UKA
ADY RÓWNAC
LINIOWYCH
INSTRUKCJA LABORATORYJNA DO ĆWICZENIA NR 4.
WERSJA ROBOCZA
OPRACOWAA
:
mgr inż. Tomasz KWAŚNIEWSKI
POLITECHNIKA ÅšWI
TOKRZYSKA
KIELCE 2013
1. Wprowadzenie
Jednym z zadań algebry liniowej jest rozwiązanie układu równań liniowych [1, 3]
= (1)
gdzie:
A  jest macierzÄ… m × n okreÅ›lonych współczynników,
 wektorem m danych liczb,
 wektorem n niewiadomych.
Układ równań opisany zależnością (1) może mieć nieskończenie wiele rozwiązań,
jedno rozwiązanie lub nie mieć ich wcale. Metody numerycznego rozwiązywania układów
równań liniowych możemy podzielić na dwie grupy [1, 3].
" metody bezpośrednie przy braku błędów zaokrągleń dają dokładne
rozwiązanie po skończonej liczbie przekształceń układu wyjściowego,
" metody iteracyjne natomiast pozwalają na wyznaczenie zbieżnego ciągu
rozwiązań przybliżonych.
Metody bezpośrednie są najbardziej efektywne w rozwiązywaniu układów o
macierzach pełnych. Ich wadami są duże obciążenie pamięci komputera oraz niestabilność
wynikająca z błędów zaokrągleń. Natomiast metody iteracyjne stosujemy układów równań
liniowych, których macierze są rzadkie.
Układ równań liniowych  postać ogólna [2]:
+ + ï" + =
+ + ï" + =
î" (2)
+ + ï" + =
Wartości m, n są dowolne. Definiujemy macierz rozszerzoną układu równań liniowych w postaci:
ï"
ï"
= (3)
î" î" î" î"
Å„"
ï"
Twierdzenie. (Capellego)
Warunkiem koniecznym i wystarczającym rozwiązywalności dowolnego układu równań liniowych jest,
aby rząd r macierzy A układu był równy rzędowi macierzy rozszerzonej D. Jeśli powyższy warunek jest
spełniony, to układ (2) posiada rozwiązanie zależne od n - r parametrów, gdzie n to liczba
niewiadomych. W przypadku n = r istnieje jednoznaczne rozwiÄ…zanie.
2. Metody bezpośrednie (dokładne)
a. Metoda Cramera
Rozwiązanie układu (2) w przypadku nieosobliwej macierzy kwadratowej A można wyznaczyć
ze wzorów Cramera [2, 3],
= , k = 1,2, & , n (4)
gdzie macierz Ak powstaje z macierzy A przez zastÄ…pienie k tej kolumny przez wektor b.
" skrypt wykorzystujÄ…cy metodÄ™ Cramera  m-plik
function x = m_cramer(A, b)
[m n] = size(b);
z = zeros(m, 1);
Ai = A;
for k=1:m
Ai(:, k) = b;
z(k) = det(Ai);
Ai(:, k) = A(:, k);
end
detA = det(A);
x = z / detA;
return
end
b. Metoda eliminacja Gaussa
Jeżeli macierz układu (1) jest macierzą trójkątną, rozwiązanie układu równań jest proste.
Przyjmijmy, że A jest macierzą trójkątną górną. Aby istniało jednoznaczne rozwiązanie, macierz musi
być nieosobliwa. Układ równań przedstawiono poniżej [1, 2, 3]:
+ + ï" + =
+ ï" + =
î" (5)
=
Składową xn wektora niewiadomych obliczymy jako pierwszą z ostatniego równania powyższego
układu. Podstawiając wynik do przedostatniego równania wyznaczymy xn-1. Procedurę kontynuujemy
aż do wyliczenia x1 . Rozwiązanie przedstawia się wzorem:
 =
(6)
"
= , = - 1, - 2, & ,1
Ponieważ obliczenia zaczynamy od ostatniej składowej wektora niewiadomych, metoda ta nazywana
jest podstawianiem w tył. W podobny sposób znajdziemy rozwiązanie układu z macierzą trójkątną
dolną, tym razem wykonując podstawianie w przód:
=
(7)
"
= , = 2,3, & ,
Wiele metod numerycznego rozwiązywania układów równań z dowolnymi macierzami polega na
sprowadzeniu układu wyjściowego do postaci trójkątnej, a następnie zastosowaniu jednego ze
wzorów (6)-(7). Jedną z takich właśnie metod jest eliminacja Gaussa. Dla ułatwienia założymy, że
macierz ukÅ‚adu jest wymiaru 3 × 3, tzn.
+ + =
+ + = (8)
+ + =
Celem jest sprowadzenie tego układu do postaci trójkątnej. Odejmując od drugiego wiersza układu
(8) pierwszy pomnożony przez a21/a11 , a od trzeciego pierwszy pomnożony przez a31/a11 otrzymamy
+ + =
+ = (9)
+ =
gdzie
= = , , = 1,2,3 (10)
,
oraz
= - , = - , , = 2,3 (11)
Aby wyznaczyć zmienną x2 z trzeciego równania w układzie (9), odejmujemy od niego drugie
pomnożone przez / :
+ + =
+ = (12)
=
gdzie
= - , = - , , = 3 (13)
Wzory na współczynniki macierzy i wyrazy wolne w każdym kroku eliminacji Gaussa łatwo jest
uogólnić na przypadek macierzy dowolnego rozmiaru:
= - , , = + 1, + 2, & , (14)
= - , , = + 1, + 2, & , (15)
Otrzymaliśmy układ trójkątny, który można teraz rozwiązać podstawianiem wstecz (6),
"
= , = , - 1, & ,1 (16)
" skrypt wykorzystujÄ…cy metodÄ™ Gaussa  m-plik [2]
function [x, det] = gauss(A, b)
if size(b, 2) > 1;
b = b';
end
n = length(b);
for k = 1:n-1
for i= k+1:n
if A(i, k) ~= 0
lambda = A(i, k)/A(k, k);
A(i, k+1:n) = A(i, k+1:n) - lambda*A(k, k+1:n);
b(i)= b(i) - lambda*b(k);
end
end
end
if nargout == 2;
det = prod(diag(A));
end
for k = n:-1:1
b(k) = (b(k) - A(k, k+1:n)*b(k+1:n))/A(k, k);
end
x = b;
c. Metoda eliminacja Jordana
Rozważmy układ równań liniowych [1, 2, 3]
+ + ï" + =
+ + ï" + =
î" (17)
+ + ï" + =
Dzielimy pierwsze równanie obustronnie przez ,a następnie od i tego wiersza (i = 2, 3, . . . , n)
odejmujemy pierwszy pomnożony przez ,
+ + ï" + =
+ ï" + =
î" (18)
+ ï" + =
W kolejnym kroku dzielimy drugie równanie obustronnie przez i odejmujemy od i tego wiersza
(i= 1, 3, 4, . . . , n) wiersz drugi pomnożony przez ,
+ ï" + =
+ ï" + =
î" (19)
ï" + =
Po (n - 1) eliminacjach otrzymujemy układ
=
=
Å„" î" (20)
=
czyli gotowe rozwiÄ…zanie.
3. Rozkład LU
Wiemy już, że rozwiązanie układu równań z macierzami trójkątnymi jest proste.
Załóżmy, że macierz A dowolnego układu da się przedstawić w postaci iloczynu macierzy
trójkątnej dolnej L i trójkątnej górnej U [1, 2, 3, 4],
A = LU (21)
Jeżeli macierz A jest nieosobliwa, zachodzi
= = (22)
a więc rozwiązanie układu (1) da się przedstawić w postaci:
= = (23)
Aby znalez
ć rozwiązanie układu dysponując rozkładem LU jego macierzy, wystarczy
rozwiązać dwa układy trójkątne:
= (24)
= (25)
a. Metoda Doolittle a
Rozkładu LU możemy poszukać również w inny sposób, traktując równość (21) [1, 2, 3, 4],
1 0 0
0 (26)
1 0
=
0 0
1
jako układ n2 równań dla n2 niewiadomych lij i uij :
+ (27)
= +
+ + +
StÄ…d
= = =
= = - = -
(28)
= = = - -
W przypadku ogólnym elementy macierzy L i U obliczamy kolejno dla i = 1, 2, & , n ze
wzorów
= - " , = , + 1, & , (29)
"
= , = + 1, + 2, & , (30)
Metoda Doolittle a staje się niezawodna dopiero w połączeniu z wyborem elementu
podstawowego. Wiersze powinniśmy zamieniać ze sobą miejscami tak, aby element uii we
wzorze (30) był jak największy.
" skrypt  dekompozycja macierzy A  wzór (24)  m-plik [2, 4]
function A = LU_Dek(A)
n = size(A, 1);
for k = 1:n-1
for i = k+1:n
if A(i,k) ~= 0.0
lambda = A(i, k)/A(k, k);
A(i, k+1:n) = A(i, k+1:n) - lambda*A(k, k+1:n);
A(i, k) = lambda;
end
end
end
" skrypt  rozwiązanie układu równań  wzór (25)  m-plik [2]
function x = LU_Solve(A, b)
if size(b,2) > 1;
b = b';
end
n = length(b);
for k = 2:n
b(k) = b(k) - A(k, 1:k-1)*b(1:k-1);
end
for k = n:-1:1
b(k) = (b(k) - A(k, k+1:n)*b(k+1:n))/A(k,k);
end
x = b;
" skrypt  wywołanie  m-plik [2]
A = [3 -1 4; -2 0 5; 7 2 -2];
B = [6 -4; 3 2; 7 -5];
A = LU_Dek(A);
det = prod(diag(A))
for i = 1:size(B,2)
X(:, i) = LU_Solve(A, B(:, i));
end
X
4. Zadania do wykonania
Zadanie 1
Znajdz
macierz A i |A| z macierzy LU
1 0 0
= 1 1 0 , = 1 2 4
0 3 21
5
1 1 0 0 0
3
2 0 0 2 -1 1
= -1 1 0 0 1 -3
, =
1 -3 1 0 0 1
Zadanie 2
MajÄ…c macierze L i U
1 0 0 2 -3 -1
3/2 1 0 0 13/2 -7/2
= = rozwiąż = gdzie = -1 2
.
1
1/2 11/13 1 0 0 32/13
Zadanie 3
Rozwiąż układ równań metodą eliminacji Gaussa =
2 -3 -1 3
= , = -9
3 2 -5
2 4 -1 -5
Zadanie 4
Znajdz
macierze L oraz U metodÄ… Doolittle s
4 -1 0
= = -1 4 -1
0 -1 4
Zadanie 5
Znajdz
rozwiązanie układu równań opisanych zależnością = przy pomocy metody
Doolittle s.
-3 6 -4 -3
= , =
9 -8 24 65
-12 24 -26 -42
2.34 -4.10 1.78 0.02
= -1.98 3.47 -2.22
, = -0.73
2.36 -15.17 6.18 -6.63
4 -3 6 1 0
= -3 10 0 1
, =
8
-4 12 -10 0 0
5. Literatura
[1] A. Jastriebow, M. Wciślik, Wstęp do metod numerycznych, Politechnika Świętokrzyska, Kielce
2000.
[2] Jaan Kiusalaas, Numerical Methods in Engineering with MATLAB, Cambridge University Press,
Cambridge 2005.
[3] Janusz Szwabiński, Metody numeryczne I, 2006.
[4] W. Y. Yang, W. Cao, T.-S. Chung, and J. Morris. Applied numerical methods using Matlab.
Wiley-Interscience, 2005.