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).
Wyszukiwarka
Podobne podstrony:
Lab 02 2011 2012
Lab 06 2011 2012
Lab 09 2011 2012
Lab 06 2011 2012 NWD
Lab 10 2011 2012
Lab 05 2011 2012
Lab 04 2011 2012
Lab 09 2011 2012
Lab 03 2011 2012
Lab 03 2011 2012
Lab 07 2011 2012 Suplement
Lab 02 2011 2012
Lab 06 2011 2012
Lab 09 2011 2012
Lab 06 2011 2012 NWD
więcej podobnych podstron