Informatyka I – Lab 08, r.a. 2011/2012
prow. Sławomir Czarnecki
Zadania na laboratorium nr. 8
Po utworzeniu nowego projektu, dołącz bibliotekę
bibs.h
.
1. Zdefiniuj funkcję
(
)
,
,
,
abx
n
void
double **
double *
double *
b
x
int
a
implementującą
Algorytm Gaussa znajdowania rozwiązania
n
∈
x
ℝ liniowego układu równań
=
a x
b
(1)
gdzie
n n
M
×
∈
a
jest daną, nieosobliwą macierzą rzędu n (
rank
n
=
a
),
n
∈
b
ℝ
jest danym
wektorem prawej strony.
Algorytm Gaussa – szkic metody
Zakładamy, Ŝe
0,0
0
a
≠
, modyfikujemy wiersze macierzy
n n
M
×
∈
a
oraz składowe wektora
n
∈
b
ℝ
o numerach
1, 2,...,
1
k
n
=
−
obliczając
,0
,0
1
1
,
,
0,
0
0,0
0,0
,
k
k
k j
k j
j
k
k
a
a
a
a
a
b
b
b
a
a
=
−
=
−
dla
1, ,...,
1
j
k
k
n
= −
−
, tzn. macierz a i wektor b
1,1
1,
1,
1
1
1,1
1,
1
1
1
0,0
0,1
0,
1
0
,0
0
...
...
,
...
...
...
...
...
...
n
n
n
n
n
n
n
a
a
b
a
a
b
a
a
a
b
a
a
−
−
−
−
−
−
−
(2)
przekształcamy do macierzy
1
a
i wektora
1
b
równych odpowiednio
0,0
0,1
0,
1
0
0,0
0,1
0,
1
0
0,0
0,0
0,0
0
1,0
1,0
1,0
1,
1,0
1,0
1,0
1,0
1,
,0
0,0
0,1
0,
1
0,
1
1,
1
1
1,1
0
0,
1
1,
0
,
0
1
0
,0
1
0
...
...
,
...
...
...
...
...
...
n
n
n
n
n
n
n
n
n
n
n
n
a
a
a
b
a
a
a
b
a
a
a
a
a
a
a
a
a
a
a
a
a
b
a
a
a
a
a
a
a
a
a
b
a
−
−
−
−
−
−
−
−
−
−
−
−
−
−
−
−
−
−
−
0
0,0
1,0
n
a
b
a
−
−
.
(3)
Otrzymany w ten sposób, zmodyfikowany układ równań
1
1
=
a x
b
(4)
charakteryzuje macierz
1
n n
M
×
∈
a
, w której wszystkie składowe w pierwszej kolumnie
poniŜej elementu
0,0
a
są równe 0:
0,0
0,1
0,
1
0
1
1
1
1,1
1,
1
1
1
1
1
1
1
1,1
1,
1
1
...
0
...
,
...
...
...
...
...
0
...
n
n
n
n
n
n
a
a
a
b
a
a
b
a
a
b
−
−
−
−
−
−
=
=
a
b
(5)
Opisane powyŜej przekształcenia, powtarzamy dla podmacierzy i wektora przedstawionych
poniŜej:
1
1
1
2,2
2,
1
2
1
1
1
1
1,1
1,2
1,
1
1
2
1
1
1
1
1,
1
1,1
2
1,
1
1
,1
...
...
...
...
...
...
...
...
...
...
,
...
...
...
...
...
...
...
...
n
n
n
n
n
n
n
a
a
b
a
a
b
a
a
a
a
b
a
−
−
−
−
−
−
−
.
(6)
Zakładając zatem, Ŝe składowa
1
1,1
0
a
≠ , modyfikujemy wiersze macierzy
1
n n
M
×
∈
a
oraz
składowe
wektora
1
n
∈
b
ℝ
o
numerach
2,3,...,
1
k
n
=
−
obliczając
1
1
,1
,1
2
1
1
2
1
1
,
,
1,
1
1
1
1,1
1,1
,
k
k
k j
k j
j
k
k
a
a
a
a
a
b
b
b
a
a
=
−
=
−
dla
1, ,...,
1
j
k
k
n
= −
−
1
1
1
1,1
1,2
1,
1
1
1
1
1,1
1,2
1,
1
1
1
1
1,1
1,1
1,1
1
1
1
1,1
1,2
1,
1
1
1
1
1,1
1,1
1,
1
1
2,2
2,
1
1
1
1,2
1
1
1
2,1
2,1
1
1
1
1,1
2,1
1
2
1,1
1,
,1
1
1,
1
1
1
1,1
...
...
...
...
...
...
...
n
n
n
n
n
n
n
n
n
n
n
a
a
a
a
a
a
a
a
a
a
a
a
a
a
a
a
a
a
a
a
a
a
a
a
a
a
a
−
−
−
−
−
−
−
−
−
−
−
−
−
−
−
−
−
1
1
2
1
1
1
1
1
1
1,1
1
1
1
1,1
1
1,
2,1
1
1
,
...
n
n
a
b
b
a
a
b
b
a
b
−
−
−
−
.
(7)
Otrzymany w ten sposób, zmodyfikowany układ równań
2
2
=
a x
b
(8)
charakteryzuje macierz
2
n n
M
×
∈
a
, w której dodatkowo, oprócz wszystkich składowych w
pierwszej kolumnie poniŜej elementu
0,0
a
, wszystkie składowe w drugiej kolumnie poniŜej
elementu
1
1,1
a są takŜe równe 0
0,0
0,1
0,2
0,
1
0
1
1
1
1
1,1
1,2
1,
1
1
2
2
2
2
2
2,2
2,
1
2
2
2
2
1,2
1,
1
1
...
0
...
0
0
...
,
...
...
...
...
...
...
0
0
...
n
n
n
n
n
n
n
a
a
a
a
b
a
a
a
b
a
a
b
a
a
b
−
−
−
−
−
−
−
=
=
a
b
(9)
Kontynuacja opisanych powyŜej obliczeń prowadzi do modyfikacji układu równań z
początkowej postaci
=
a x
b
do postaci
1
1
n
n
−
−
=
a
x
b
(10)
(z tzw. trójkątną macierzą układu równań), gdzie
0,0
0,1
0,2
0,
1
0
1
1
1
1
1,1
1,2
1,
1
1
2
2
1
1
2
2,2
2,
1
2
1
1
1,
1
1
...
0
...
0
0
...
,
...
...
...
...
...
...
0
0
0
...
n
n
n
n
n
n
n
n
n
n
a
a
a
a
b
a
a
a
b
a
a
b
a
b
−
−
−
−
−
−
−
−
−
−
=
=
a
b
.
(11)
W przypadku gdyby któraś ze składowych
0,0
a
,
1
1,1
a ,
2
2,2
a
,… była równa 0, konieczna jest
dodatkowa zmiana lub zmiany polegające na odszukaniu (w najprostszym wariancie)
pierwszego wiersza k („poniŜej wiersza
1
i −
”), w którym składowa
(
)
,0
0
k
a
k >
,
(
)
1
,1
1
k
a
k >
,
(
)
2
,2
2
k
a
k >
,… jest róŜna od 0 i zamianie wiersza, w którym jest zerowa składowa ze
znalezionym wierszem o numerze k. W praktyce, znajdowanie odpowiednich składowych jest
nieco bardziej złoŜone i wynika z konieczności zapewnienia minimalizacji błędów
numerycznych, które sukcesywnie pojawiają się w trakcie przeprowadzania opisanych
powyŜej obliczeń.
Oznaczając na nowo macierz
1
n
n n
M
−
×
∈
a
jaki i wszystkie jej składowe
1
,
n
i j
a
−
(
)
,
0,1,...,
1
i j
n
=
− , a takŜe wektor
1
n
n
−
∈
b
ℝ
jak i wszystkie jego składowe
1
n
i
b
−
(
)
0,1,...,
1
i
n
=
− ponownie symbolami odpowiednio
( )
,
i j
n n
a
M
×
=
∈
a
oraz
( )
n
i
b
=
∈
b
ℝ ,
stwierdzamy, Ŝe zadanie poszukiwania rozwiązania układu równań liniowych (1)
sprowadzone zostało do zadania
=
a x
b
(12)
z trójkątną macierzą
n n
M
×
∈
a
, gdzie
0,0
0,1
0,
1
1,1
1,
1
1,
1
...
0
...
...
...
...
...
0
0
...
n
n
n n
n
n
a
a
a
a
a
M
a
−
−
×
−
−
=
∈
a
,
0
1
1
...
n
n
b
b
b
−
=
∈
b
ℝ .
(13)
Rozwiązanie
( )
n
i
x
=
∈
x
ℝ układu równań (12) moŜna w bardzo prosty sposób znaleźć
wykorzystując trójkątną postać macierzy a. Najpierw znajdujemy ostatnią składową
1
n
x
−
z
ostatniego równania (13),
1
1
1,
1
n
n
n
n
b
x
a
−
−
−
−
=
, a kolejne składowe obliczamy rekurencyjnie ze
wzoru
(
)
1
,
1
,
2,
3,...,1, 0
n
i
i j
j
j i
i
i i
b
a
x
x
i
n
n
a
−
= +
−
=
= −
−
∑
.
Algorytm Gaussa (najprostsza wersja) – pseudo-kod
start
dla kaŜdego wiersza
1, 2,...,
1
i
n
=
−
{
jeśli
1, 1
0
i
i
a
− −
=
, to znajdź pierwszy wiersz
(
)
1
k i
k
n
≤ ≤ − , w którym
, 1
0
k i
a
−
≠
, a
następnie zamień wiersz
1
i − z wierszem k aktualnej macierzy a oraz
zamień składową
1
i − ze składową k aktualnego wektora b, tzn.
(
)
1,
,
1
1, ,...,
1 ,
i
j
k j
i
k
a
a
j
i
i
n
b
b
−
−
↔
= −
−
↔
(moŜna ograniczyć się do zamiany składowych macierzy a jedynie od
1
i − kolumny)
dla kaŜdego wiersza
(
)
1
k i
k
n
≤ ≤ −
{
, 1
1, 1
k i
i
i
a
c
a
−
− −
=
, 1
0
k i
a
−
=
dla kaŜdej kolumny
(
)
1
j i
j
n
≤ ≤ −
{
,
,
1,
k j
k j
i
j
a
a
c a
−
=
−
}
1
k
k
i
b
b
c b
−
=
−
}
}
1
1
1,
1
n
n
n
n
b
x
a
−
−
−
−
=
dla kaŜdej składowej
2,
3,...,1, 0
i
n
n
= −
−
{
1
,
1
,
n
i
i j
j
j i
i
i i
b
a
x
x
a
−
= +
−
=
∑
}
koniec
2. Przetestuj działanie funkcji
(
)
,
,
,
abx
n
void
double **
double *
double *
b
x
int
a
dla
macierzy
N N
M
×
∈
A
i wektora
N
∈
B
ℝ o losowo wygenerowanych składowych z przedziału
[
]
10,10
−
oraz
10
N =
. Wyświetl na ekranie rozwiązanie
N
∈
X
ℝ układu równań
=
A X
B
, a
takŜe wygenerowaną macierz
N N
M
×
∈
A
oraz wektor
N
∈
B
ℝ przed i po wywołaniu funkcji
( )
...
abx
void
.
3. Dla wygenerowanej w zadaniu 2 macierzy
N N
M
×
∈
A
i wektora
N
∈
B
ℝ , wywołaj funkcję
(
)
,
,
,
gauss
n
int
int
double **
a
double *
d
b
ouble *
x z biblioteki
bibs.h
. Porównaj wynik z
rozwiązaniem numerycznym
N
∈
X
ℝ z zadania 2.
4. Przetestuj działanie funkcji
(
)
,
,
,
abx
n
void
double **
double *
double *
b
x
int
a
dla
następująco zdefiniowanej macierzy
2 2
M
×
∈
A
i wektora
2
∈
B
ℝ (N = 2)
2
2 2
1
1
,
1 1
2
M
ε
×
=
∈
=
∈
A
B
ℝ ,
(14)
gdzie
ε
jest „bardzo małą” liczbą dodatnią, np.:
20
1.0 10
ε
−
=
⋅
.
Komentarz. Rozwiązanie ścisłe (dokładne)
0
2
1
X
X
=
∈
X
ℝ układu równań
=
A X
B
(15)
tzn. układu
0
1
1
1
1 1
2
X
X
ε
=
(16)
jest następujące
1
0
1
1
2
1
1
1
1
...
...
1
1
1 2
1
2
1
1
1
1
1
1
2
X
X
X
ε
ε
ε
ε
ε
ε
ε
ε
ε
ε
ε
−
−
−
−
−
−
=
=
= =
= =
−
−
−
−
−
X
.
(17)
Algorytm Gaussa (w najprostszej postaci) doprowadza układ równań (16) do układu równań
następującego:
0
1
1
1
1
1
0 1
2
X
X
ε
ε
ε
=
−
−
,
(18)
który ma rozwiązanie ścisłe (dokładne)
1
1
1
2
1
X
ε
ε
=
−
−
,
1
0
1 X
X
ε
−
=
.
(19)
WaŜne !
Kolejność
1
0
,
X
X zapisanych w (19) składowych odpowiada kolejności w jakiej
obliczane są one w algorytmie Gaussa – patrz pseudo-kod algorytmu Gaussa wyŜej.
Dla bardzo małego
0
ε
> , np.
20
1.0 10
ε
−
=
⋅
, wartości obu wyraŜeń
1
2
ε
−
oraz
1
1
ε
−
w (19)
1
obliczanych na komputerze są w przybliŜeniu równe
1
ε
−
, tzn. zarówno
2
1
1
ε
ε
−
≈
−
jak i
1
1
1
ε
ε
−
≈
−
. A zatem ze wzoru (19)
1
otrzymamy następujący wynik numeryczny
1
1
1
2
1
1
1
1
X
ε
ε
ε
ε
−
−
−
−
=
≈
= ,
(20)
i w konsekwencji ze wzoru (19)
2
obliczymy
1
0
1
1 1
0
X
X
ε
ε
−
−
=
≈
=
.
(21)
Jednak, po prostym przekształceniu wzorów (19) (porównaj takŜe (17)) otrzymujemy
następujące, poprawne numerycznie wartości:
1
2
1
1
1
X
ε
ε
−
−
=
≈
,
0
1
1
1
X
ε
=
−
≈
.
(22)
Otrzymane rozwiązanie numeryczne otrzymane najprostszą metodą Gaussa jest zatem nie do
zaakceptowania. Przyczyną, tak duŜego, numerycznego błędu jest fakt, Ŝe składowa
0,0
a
ε
=
jest „bardzo mała” w porównaniu z innymi składowymi pierwszego wiersza (u nas w
porównaniu ze składową
0,1
1
a
=
), co nie moŜe być w prosty sposób poprawnie
zaimplementowane na komputerach o stałej długości słowa, w którym musi być
przechowywana reprezentacja liczb zmienno przecinkowych w postaci zarówno cechy jak i
mantysy.
ZauwaŜmy jednak, Ŝe proste przestawienie kolejności równań (16)
0
1
1 1
2
1
1
X
X
ε
=
(23)
a następnie zastosowanie algorytmu eliminacji Gaussa (w najprostszej wersji) prowadzi do
równowaŜnego układu równań
0
1
1
1
2
0 1
1 2
X
X
ε
ε
=
−
−
,
(24)
który tym razem, ma poprawne numerycznie rozwiązanie
1
1 2
1
1
X
ε
ε
−
=
≈
−
(25)
oraz
0
1
2
1
X
x
= −
≈ .
(26)
Wynika stąd, Ŝe dobry pod względem numerycznym algorytm rozwiązywania układu
równań liniowych powinien w ogólnym przypadku brać pod uwagę moŜliwość m.in.
odpowiedniego przestawienia wierszy w celu eliminacji przedstawionych wyŜej
problemów numerycznych.
5. Raz jeszcze przetestuj działanie funkcji
( )
...
abx
void
, ale tym razem dla następująco
zdefiniowanej macierzy
2 2
M
×
∈
A
i wektora
2
∈
B
ℝ (N = 2)
2
2 2
1 1
2
,
1
1
M
ε
×
=
∈
=
∈
A
B
ℝ
(27)
(prosta zamiana kolejności równań (16)). Wyświetl na ekranie rozwiązanie
N
∈
X
ℝ .
Dla zdefiniowanej w (27) macierzy
N N
M
×
∈
A
i wektora
N
∈
B
ℝ , wywołaj raz jeszcze
funkcję
(
)
,
,
,
gauss
n
int
int
double **
a
double *
d
b
ouble *
x z biblioteki
bibs.h
i wyświetl
na ekranie rozwiązanie
N
∈
X
ℝ . Porównaj oba otrzymane wyniki z (niepoprawnym)
rozwiązaniem numerycznym (20), (21) oraz z (poprawnym) rozwiązaniem (25), (26).