Rozwi ˛
azywanie układów równa´n liniowych
1
4 Rozwi ˛
azywanie układów równa ´n liniowych
Standardowa posta´c układu równa ´n:
A
1,1
A
1,2
· · ·
A
1,n−1
A
n,n
A
2,1
A
2,2
· · ·
A
2,n−1
A
n,n
..
.
..
.
..
.
..
.
..
.
A
n−1,1
A
n−1,2
· · · A
n−1,n−1
A
n−,n
A
n,1
A
n,2
· · ·
A
n,n−1
A
n,n
|
{z
}
A
·
x
1
x
2
..
.
x
n−1
x
n
|
{z
}
x
=
b
1
b
2
..
.
b
m−1
b
m
|
{z
}
b
Z algebry wiadomo, ˙ze je´sli macierz
A
jest
nieosobliwa
, to jedyne rozwi ˛
azanie tego układu
wyra˙za si ˛e jako:
x = A
−1
· b
, gdzie
A
−1
jest macierz ˛
a odwrotn ˛
a do macierzy
A
.
L.J. Opalski, “Wst ˛ep do metod numerycznych - materiały do wykładu”
wersja: 16.III.2010
Rozwi ˛
azywanie układów równa´n liniowych
2
Wyznaczanie rozwi ˛
azania układu równa ´n
x = A
−1
· b
w dwóch krokach:
1. oblicz macierz
C
odwrotn ˛
a do
A
, t.j.:
C = A
−1
2. oblicz rozwi ˛
azanie:
x = C · b
nie jest korzystne
ze wzgl ˛edu na
nakłady obliczeniowe i dokładno´s´c
.
Uzasadnienie-przykład
Równanie
7 · x = 21
ma rozwi ˛
azanie całkowite:
x = 21/7 = 3
. Wyznaczanie rozwi ˛
azania
przy u˙zyciu odwrotno´sci współczynnika przy
x
, tzn.:
x = 7
−1
· 21 ≈ 0.142857 · 21 = 2.99997
, jest mniej dokładne i wymaga wi ˛ecej operacji
arytmetycznych.
Dalej b ˛edzie mowa:
•
o nieiteracyjnych algorytmach efektywnego i dokładnego rozwi ˛
azywania równa ´n liniowych
•
o tym jak na
dokładno´s´c wyznaczania rozwi ˛
azania
wpływa
dokładno´s´c danych
(elementów
macierzy
A
i wektora
b
) oraz
wybór algorytmu
.
•
o algorytmach iteracyjnych i obszarach ich zastosowa ´n
L.J. Opalski, “Wst ˛ep do metod numerycznych - materiały do wykładu”
wersja: 16.III.2010
Rozwi ˛
azywanie układów równa´n liniowych
3
4.1 Nieiteracyjne rozwi ˛
azywanie układów równa ´n
4.1.1 Metoda eliminacji Gaussa
Eliminacja Gaussa jest algorytmem dwufazowym:
•
Pierwsza faza przekształca obydwie strony układu równa ´n tak, by uzyska´c
równowa˙zny
układ z macierz ˛
a górn ˛
a trójk ˛
atn ˛
a (
U
):
U · x = z
W ka˙zdym kroku przekształce ´n zerowane s ˛
a elementy kolejnej kolumny macierzy układu
poni˙zej diagonali
•
W fazie drugiej rozwi ˛
azywany jest układ równa ´n
U · x = z
wzgl ˛edem kolejnych
zmiennych - pocz ˛
awszy od ostatniej.
L.J. Opalski, “Wst ˛ep do metod numerycznych - materiały do wykładu”
wersja: 16.III.2010
Rozwi ˛
azywanie układów równa´n liniowych
4
Przykład
R
1
:
R
2
:
R
3
:
3
1
6
2
1
3
1
1
1
|
{z
}
A
x
1
x
2
x
3
| {z }
x
=
2
7
4
| {z }
b
=⇒
R
1
bez zmian
R
2
:= R
2
− R
1
∗
2
/
3
R
3
:= R
3
− R
1
∗
1
/
3
=⇒
1
0
0
2/3
1
0
1/3
0
1
|
{z
}
L
(1)
R
1
:
R
2
:
R
3
:
3
1
6
0
1
3
−1
0
2
3
−1
|
{z
}
A
(1)
x
1
x
2
x
3
| {z }
x
=
2
17
3
10
3
| {z }
b
(1)
=⇒
R
1
bez zmian
R
2
bez zmian
R
3
:= R
3
− R
2
∗
2
3
1
3
=⇒
1
0
0
2
3
1
0
1
3
2
1
|
{z
}
L≡L
(2)
R
1
:
R
2
:
R
3
:
3
1
6
0
1
3
-1
0
0
1
|
{z
}
U≡A
(2)
x
1
x
2
x
3
| {z }
x
=
2
17
3
−8
| {z }
z≡b
(2)
Uwaga:
L
(1)
A
(1)
= A
,
L
(1)
b
(1)
= b
L
(2)
A
(2)
= A
,
L
(2)
b
(2)
= b
LU = A
,
Lz = b
R
1
:
x
1
=
2−
1
·x
2
−
6
·x
3
3
R
2
:
x
2
=
17/3−(
-1
)·x
3
1/3
R
3
:
x
3
=
−8
1
⇒
x =
19
−7
−8
L.J. Opalski, “Wst ˛ep do metod numerycznych - materiały do wykładu”
wersja: 16.III.2010
Rozwi ˛
azywanie układów równa´n liniowych
5
Przykład
Znajd´zmy napi ˛ecia w˛ezłowe
U
1
,
U
2
,
U
3
,
U
4
w elektrycznej sieci liniowej:
G
34
G
3
G
23
G
12
G
2
G
1
¹¸
º·
±°
²¯
6
I = 1
1
2
3
4
Równania bilansu pr ˛
adów (dla w˛ezłow 1, 2, 3, 4)
daj ˛
a:
G
1
U
1
+ G
12
(U
1
− U
2
)
=
0
G
2
U
2
+ G
12
(U
2
− U
1
) + G
23
(U
2
− U
3
)
=
0
G
3
U
3
+ G
23
(U
3
− U
2
) + G
3
(U
3
− U
4
)
=
0
G
4
(U
4
− U
3
)
=
1
Przyjmijmy
G
i
= 1
,
G
ij
= 2
, otrzymuj ˛
ac
układ równa ´n:
3
−2
0
0
−2
5
−2
0
0
−2
5
−2
0
0
−2
2
U
1
U
2
U
3
U
4
=
0
0
0
1
Sie´c elektryczn ˛
a
2
1
2
2
1
1
¹¸
º·
±°
²¯
6
1
1
2
3
4
-
G
eq,1
=
2·1
2+1
=
2
3
L.J. Opalski, “Wst ˛ep do metod numerycznych - materiały do wykładu”
wersja: 16.III.2010
Rozwi ˛
azywanie układów równa´n liniowych
6
Etap I.
Eliminacja zmiennych
1. Mno˙z ˛
ac równanie 1 przez 2/3 i dodaj ˛
ac do drugiego:
3
−2
0
0
0
11/3
−2
0
0
−2
5
−2
0
0
−2
2
U
1
U
2
U
3
U
4
=
0
0
0
1
2
1
2
5
3
¹¸
º·
±°
²¯
6
1
2
3
4
-
G
eq,2
=
2·5/3
2+5/3
=
10
11
2. Mno˙z ˛
ac równanie 2 przez 6/11 i dodaj ˛
ac do trzeciego:
3
−2
0
0
0
11/3
−2
0
0
0
43/11
−2
0
0
−2
2
U
1
U
2
U
3
U
4
=
0
0
0
1
2
21
11
¹¸
º·
±°
²¯
6
1
3
4
-
G
eq,3
=
2·21/11
2+21/11
=
42
43
L.J. Opalski, “Wst ˛ep do metod numerycznych - materiały do wykładu”
wersja: 16.III.2010
Rozwi ˛
azywanie układów równa´n liniowych
7
3. Mno˙z ˛
ac równanie 3 przez 22/43 i dodaj ˛
ac do czwartego:
3
−2
0
0
0
11/3
−2
0
0
0
43/11
−2
0
0
0
42/43
U
1
U
2
U
3
U
4
=
0
0
0
1
42
43
¹¸
º·
±°
²¯
6
1
4
Etap II.
Podstawienia odwrotne
(wstecz)
3
−2
0
0
0
11/3
−2
0
0
0
43/11
−2
0
0
0
42/43
U
1
U
2
U
3
U
4
=
0
0
0
1
St ˛
ad:
42
43
U
4
=
1 ⇒
U
4
=
43
42
43
11
U
3
− 2U
4
= 0
⇒
U
3
=
11
21
11
3
U
2
− 2U
3
= 0
⇒
U
2
=
6
21
3U
1
− 2U
2
= 0
⇒
U
1
=
4
21
L.J. Opalski, “Wst ˛ep do metod numerycznych - materiały do wykładu”
wersja: 16.III.2010
Rozwi ˛
azywanie układów równa´n liniowych
8
Wykonalno ´s ´c algorytmu Gaussa
Przykład
Układ równa ´n
1 2 3
1 2 1
2 1 3
|
{z
}
A
x =
6
4
6
| {z }
b
ma dokładnie jedno rozwi ˛
azanie:
1
1
1
Tymczasem po jednym kroku eliminacji Gaussa mamy:
1
2
3
0
0
−2
0 −3 −3
x =
6
−2
−6
Drugiego kroku eliminacji
nie mo˙zna zrobi´c
!
L.J. Opalski, “Wst ˛ep do metod numerycznych - materiały do wykładu”
wersja: 16.III.2010
Rozwi ˛
azywanie układów równa´n liniowych
9
Rozwi ˛
azanie problemu
: przestawienie równa ´n 2 i 3.
1 0 0
0 0 1
0 1 0
|
{z
}
P
1 2 3
1 2 1
2 1 3
|
{z
}
A
|
{z
}
A
m
x =
1 0 0
0 0 1
0 1 0
|
{z
}
P
6
4
6
| {z }
b
|
{z
}
b
m
gdzie
P
jest
macierz ˛
a przestawie ´n
.
A
m
=
1 2 3
2 1 3
1 2 1
, b
m
=
6
6
4
Po jednym kroku eliminacji Gaussa mamy układ równa ´n z macierz ˛
a górn ˛
a trójk ˛
atn ˛
a:
1
2
3
0 −3 −3
0
0
−2
|
{z
}
U
x =
6
−6
−2
| {z }
z
L.J. Opalski, “Wst ˛ep do metod numerycznych - materiały do wykładu”
wersja: 16.III.2010
Rozwi ˛
azywanie układów równa´n liniowych
10
Przestawienia równa ´n wpływaj ˛
a nie tylko na wykonalno´s´c eliminacji Gaussa, ale i na
dokładno´s´c
Przykład
Rozwi ˛
aza´c układ równa ´n korzystaj ˛
ac z aryt-
mometru wykonuj ˛
acego operacje zmiennopozy-
cyjne z 5 dziesi ˛etnymi cyframi znacz ˛
acymi.
10
−7 0
−3 2.099 6
5
−1 5
x
1
x
2
x
3
=
7
3.901
6
Po pierwszym kroku eliminacji Gaussa:
10
−7 0
0
-0.001
6
0
2.5 5
x
1
x
2
x
3
=
7
6.001
2.5
Po drugim kroku eliminacji Gaussa:
10
−7
0
0
−0.001
6
0
0
15005
x
1
x
2
x
3
=
7
6.001
15006
| {z }
z
3
Podstawienia wstecz daj ˛
a:
x
3
= 15006/15005 = 1.0001
,
x
2
= (6.001 − 6 ∗ x
3
)/(−0.001) = −0.4
,
x
1
= (7 − 7 · x
2
)/10 = 0.9800
.
Warto´sci dokładne:
x
1
= 0
,
x
2
= −1
,
x
3
= 1
.
Przyczyna problemu: zły wybór elementu centralnego w kroku 1
⇒
du˙za wra˙zliwo´s´c na bł ˛
ad
zaokr ˛
agle ´n przy obliczaniu
z
3
.
L.J. Opalski, “Wst ˛ep do metod numerycznych - materiały do wykładu”
wersja: 16.III.2010
Rozwi ˛
azywanie układów równa´n liniowych
11
Posta´c układu równa ´n po zamianie dru-
giego i trzeciego równania
10
−7 0
5
−1 5
−3 2.099 6
x
1
x
2
x
3
=
7
6
3.901
Po pierwszym kroku eliminacji Gaussa:
10
−7 0
0
2.5 5
0 −0.001 6
x
1
x
2
x
3
=
7
2.5
6.001
Po drugim kroku eliminacji Gaussa:
10 −7
0
0 2.5
5
0
0 6.002
x
1
x
2
x
3
=
7
2.5
6.002
Podstawienia wstecz daj ˛
a
warto´sci dokładne
:
x
3
= 1
,
x
2
= (2.5 − 5 · x
3
)/2.5 = −1
,
x
1
= (7 + 7 ∗ x
2
)/10 = 0
.
L.J. Opalski, “Wst ˛ep do metod numerycznych - materiały do wykładu”
wersja: 16.III.2010
Rozwi ˛
azywanie układów równa´n liniowych
12
Uwagi
•
Przestawianie wierszy (b ˛
ad´z wierszy i kolumn), prowadz ˛
ace do du˙zych (bezwzgl ˛ednie)
warto´sci elementów centralnych (diagonalnych) zmniejsza przenoszenie bł ˛edów
zaokr ˛
agle ´n.
•
W dwóch wa˙znych przypadkach mo˙zna nie wykonywa´c przestawie ´n:
– je˙zeli macierz główna
A
jest symetryczna i dodatnio okre´slona, tzn.
A
T
= A,
i
x
T
Ax > 0,
dla ka˙zdego
x 6= 0
– je˙zeli macierz główna ma dominuj ˛
ac ˛
a przek ˛
atn ˛
a, tzn:
|a
i,i
| >
n
X
j=1,j6=i
|a
i,j
|
L.J. Opalski, “Wst ˛ep do metod numerycznych - materiały do wykładu”
wersja: 16.III.2010
Rozwi ˛
azywanie układów równa´n liniowych
13
Nakłady obliczeniowe algorytmu eliminacji Gaussa
mno˙zenie i dzielenie:
M =
1
3
n
3
+ n
2
−
1
3
n
dodawanie i odejmowanie:
D =
1
3
n
3
+
1
2
n
2
−
5
6
n
Uwagi
•
wzory Cramera wymagaj ˛
a
(n + 1)!
mno˙ze ´n
• MATLAB
: rozwi ˛
azanie układu równa ´n
Ax = b
za pomoc ˛
a:
x = A\b
L.J. Opalski, “Wst ˛ep do metod numerycznych - materiały do wykładu”
wersja: 16.III.2010
Rozwi ˛
azywanie układów równa´n liniowych
14
4.1.2 Wykorzystanie rozkładu LU
Załó˙zmy, ˙ze
A = L · U
gdzie:
L
jest
macierz ˛
a doln ˛
a trójk ˛
atn ˛
a
,
U
jest
macierz ˛
a górn ˛
a trójk ˛
atn ˛
a
.
Rozwi ˛
azanie układu równa ´n
A · x = b
uzyskuje si ˛e w dwóch krokach:
1.
podstawienia w przód
:
L · z = b
, gdzie
z
- wektor pomocniczy
2.
podstawienia wstecz
:
U · x = z
L.J. Opalski, “Wst ˛ep do metod numerycznych - materiały do wykładu”
wersja: 16.III.2010
Rozwi ˛
azywanie układów równa´n liniowych
15
Przykład
1
1
1
1 −1 1
1
2
3
|
{z
}
A
x
1
x
2
x
3
| {z }
x
=
3
1
6
| {z }
b
, A =
1
0
0
1
1
0
1 −1/2 1
|
{z
}
L
·
1
1
1
0 −2 0
0
0
2
|
{z
}
U
1
0
0
1
1
0
1 −1/2 1
|
{z
}
L
·
1
1
1
0 −2 0
0
0
2
|
{z
}
U
x
1
x
2
x
3
| {z }
x
|
{z
}
z
=
3
1
6
| {z }
b
,
1
0
0
1
1
0
1 −1/2 1
|
{z
}
L
·
z
1
z
2
z
3
| {z }
z
=
3
1
6
| {z }
b
=⇒
z
1
z
2
z
3
| {z }
z
=
3
−2
2
1
1
1
0 −2 0
0
0
2
|
{z
}
U
x
1
x
2
x
3
| {z }
x
=
3
−2
2
| {z }
z
=⇒
x
1
x
2
x
3
| {z }
x
=
1
1
1
L.J. Opalski, “Wst ˛ep do metod numerycznych - materiały do wykładu”
wersja: 16.III.2010
Rozwi ˛
azywanie układów równa´n liniowych
16
Nakłady obliczeniowe dla metody LU
rozkład macierzy:
mno˙zenie i dzielenie:
M =
1
3
n
3
−
1
3
n
dodawanie i odejmowanie:
D =
1
3
n
3
−
1
2
n
2
+
1
6
n
podstawienia:
mno˙zenie i dzielenie:
M = n
2
dodawanie i odejmowanie:
D = n
2
− n
Uwaga
Nakłady obliczeniowe na rozwi ˛
azanie układu równa ´n za pomoc ˛
a metody LU s ˛
a
takie same
jak
dla eliminacji Gaussa.
L.J. Opalski, “Wst ˛ep do metod numerycznych - materiały do wykładu”
wersja: 16.III.2010
Rozwi ˛
azywanie układów równa´n liniowych
17
Rozkład LU w programie
MATLAB
1.
[L,U,P]=lu(A)
rozkłada macierz
A
tak, ˙ze
PA = LU
, gdzie
L
jest macierz ˛
a
trójk ˛
atn ˛
a doln ˛
a,
U
jest macierz ˛
a trójk ˛
atn ˛
a górn ˛
a, a
P
jest macierz ˛
a przestawie ´n.
Ax = b
⇒
PA
|{z}
LU
x = Pb
⇒
L Ux
|{z}
z
= Pb
Wektor
x
wynika z rozwi ˛
azania dwóch układów z trójk ˛
atnymi macierzami:
Lz = Pb, Ux = z
2.
[Lm,U]=lu(A)
oblicza rozkład macierzy
A
na iloczyn
A = L
m
U
, przy czym
U
jest
macierz ˛
a trójk ˛
atn ˛
a górn ˛
a, a
L
m
jest iloczynem macierzy przestawie ´n i macierzy dolnej
trójk ˛
atnej. Rozwi ˛
azanie układu równa ´n
Ax = b
uzyskuje si ˛e z dwóch układów z
trójk ˛
atnymi macierzami:
L
m
z = b, Ux = z
L.J. Opalski, “Wst ˛ep do metod numerycznych - materiały do wykładu”
wersja: 16.III.2010
Rozwi ˛
azywanie układów równa´n liniowych
18
Przykład
1
1
1
1
2
3
1 −1 1
x
1
x
2
x
3
=
3
6
1
[ Lm,Um] = lu (A)
% Lm =
%
1.0000
0
0
%
1.0000
−
0.5000
1.0000
%
1.0000
1.0000
0
% Um =
%
1
1
1
%
0
−
2
0
%
0
0
2
z2=Lm\ b ;
x2=Um\ z2
%
1
%
1
%
1
norm ( x2
−
ones ( 3 , 1 ) , i n f )
% 0
A= [ 1 , 1 , 1 ; 1 , 2 , 3 ; 1 ,
−
1 , 1 ] ;
b = [ 3 ; 1 ; 6 ] ;
[ L , U, P] = lu (A)
% L =
%
1.0000
0
0
%
1.0000
1.0000
0
%
1.0000
−
0.5000
1.0000
% U =
%
1
1
1
%
0
−
2
0
%
0
0
2
% P =
%
1
0
0
%
0
0
1
%
0
1
0
z1=L \ ( P
∗
b ) ;
x1=U\ z1 ;
norm ( x1
−
ones ( 3 , 1 ) , i n f )
%
0
Uwaga: wyj ˛
atkowa sytuacja - brak bł ˛edów zaokr ˛
agle ´n. Powtórzy´c obliczenia, zmieniaj ˛
ac:
A
2,1
= 1.1, b
2
= 6.1
.
L.J. Opalski, “Wst ˛ep do metod numerycznych - materiały do wykładu”
wersja: 16.III.2010
Rozwi ˛
azywanie układów równa´n liniowych
19
Przykład
Znale´z´c transmitancj ˛e napi ˛eciow ˛
a
K(jω) = U
2
(jω)/E(jω)
liniowej sieci elektrycznej, jak na rysunku, dla warto´sci
parametrów
R
1
= 1kΩ, R
2
= 1M Ω, C
1
= C
2
= 1µF
. Wykre´sli´c moduł i faz˛e transmitancji dla
ω ∈ [0.1, 10
5
]
.
e(t)
C
R
R
u
u
C
1
2
1
2
1
2
Równania bilansu pr ˛
adów:
(E − U
1
)/R
1
+ (U
2
− U
1
)/R
2
=
U
1
· jωC
1
(U
1
− U
2
)/R
2
=
U
2
· jωC
2
Układ równa ´n liniowych:
1/R
1
+ 1/R
2
+ jωC
1
−1/R
2
−1/R
2
1/R
2
+ jωC
2
U
1
U
2
=
E/R
1
0
R1=1e3 ; R2=1e6 ; C1=1e
−
6; C2=C1 ;
AR= [ 1 / R1+1/R2,
−
1/R2;
−
1/R2 , 1 / R2 ] ;
AI = [C1 , 0 ; 0 , C2 ] ; b = [ 1 / R1 ; 0 ] ;
omega=logspace (
−
1 ,5);
f o r n =1: numel ( omega )
A=AR+ j
∗
omega ( n)
∗
AI ;
U=A \ b ;
T ( n )=U ( 2 ) ;
end
subplot ( 2 , 1 , 1 ) ;
semilogx ( omega,20
∗
log10 ( abs ( T ) ) ) ; axis t i g h t ; grid on ;
y l a b e l (
’|T| [dB]’
) ;
x l a b e l (
’\omega [rad/s]’
) ;
t i t l e (
’transmitancja T(j\omega)’
) ;
subplot ( 2 , 1 , 2 ) ; axis t i g h t ; grid on ;
semilogx ( omega , angle ( T ) ) ;
y l a b e l (
’arg(T) [rad]’
) ; x l a b e l (
’\omega [rad/s]’
) ;
10
0
10
2
10
4
−140
−120
−100
−80
−60
−40
−20
|T| [dB]
ω
[rad/s]
transmitancja T(j
ω
)
10
0
10
2
10
4
−3
−2
−1
arg(T) [rad]
ω
[rad/s]
L.J. Opalski, “Wst ˛ep do metod numerycznych - materiały do wykładu”
wersja: 16.III.2010
Rozwi ˛
azywanie układów równa´n liniowych
20
4.1.3 Rozkład QR
A = Q · R
R
−
macierz trójk ˛
atna górn ˛
a
Q
−
macierz ortogonalna
Poniewa˙z
Q
−1
= Q
T
, st ˛
ad:
A · x = b
Q · R · x = b
R · x = Q
T
· b
L.J. Opalski, “Wst ˛ep do metod numerycznych - materiały do wykładu”
wersja: 16.III.2010
Rozwi ˛
azywanie układów równa´n liniowych
21
Nakłady obliczeniowe metody QR
rozkład macierzy
mno˙zenie i dzielenie:
M
≈
2
3
n
3
dodawanie i odejmowanie:
D ≈
2
3
n
3
pierwiastki kwadratowe:
P
= n − 1
rozwi ˛
azywanie
mno˙zenie i dzielenie:
M
≈
3
2
n
2
dodawanie i odejmowanie:
D ≈
3
2
n
2
Rozkład QR w programie
MATLAB
1.
[Q,R]=qr(A)
rozkłada macierz
A
na iloczyn
A=QR
, przy czym
R
jest macierz ˛
a trójk ˛
atn ˛
a górn ˛
a, a
Q
macierz ˛
a unitarn ˛
a.
2.
[Q,R,P]=qr(A)
zwraca dodatkowo macierz przestawie ´n kolumn
P
tak ˛
a, ˙ze
AP=QR
, a elementy wektora
abs(diag(R))
s ˛
a malej ˛
ace.
L.J. Opalski, “Wst ˛ep do metod numerycznych - materiały do wykładu”
wersja: 16.III.2010
Rozwi ˛
azywanie układów równa´n liniowych
22
Przykład
A= [ 1 , 1 , 1 ; 1 , 2 , 3 ; 1 ,
−
1 , 1 ] ;
b = [ 3 ; 1 ; 6 ] ;
[Q,R] = qr (A)
% Q =
%
−
0.5774
−
0.1543
−
0.8018
%
−
0.5774
−
0.6172
0.5345
%
−
0.5774
0.7715
0.2673
% R =
%
−
1.7321
−
1.1547
−
2.8868
%
0
−
2.1602
−
1.2344
%
0
0
1.0690
z=Q’
∗
b ;
x=R\ z
% x =
%
1.0000
%
1.0000
%
1.0000
norm ( x
−
ones ( 3 , 1 ) , i n f )
%
1.7764e
−
015
[Q, R, P] = qr (A)
% Q =
%
−
0.3015
−
0.2752
−
0.9129
%
−
0.9045
−
0.2202
0.3651
%
−
0.3015
0.9358
−
0.1826
% R =
%
−
3.3166
−
1.8091
−
1.5076
%
0
−
1.6514
0.4404
%
0
0
−
0.7303
% P =
%
0
0
1
%
0
1
0
%
1
0
0
z=Q’
∗
b ;
x=P
∗
(R\ z )
% x =
%
1.0000
%
1.0000
%
1.0000
norm ( x2
−
ones ( 3 , 1 ) , i n f )
%
1.5543e
−
015
Uwaga: bł ˛edy zaokr ˛
agle ´n przy rozwi ˛
azywaniu równa ´n dały rozwi ˛
azanie o niedokładno´sci 14 EPS.
L.J. Opalski, “Wst ˛ep do metod numerycznych - materiały do wykładu”
wersja: 16.III.2010
Rozwi ˛
azywanie układów równa´n liniowych
23
4.2 Dokładno ´s ´c rozwi ˛
azywania układów równa ´n liniowych
4.2.1 Normy wektorów i macierzy
Wła´sciwo´sci normy wektora
k · k
w przestrzeni liniowej
R
n
:
kxk > 0 ∀x ∈ R
n
,
kxk = 0 ⇔ x = 0
kαxk = |α| · kxk, ∀α ∈ R, ∀x ∈ R
n
kx + yk 6 kxk + kyk, ∀x, y ∈ R
n
Normy Höldera (
p
-normy) wektora z liniowej przestrzeni wektorów
N
-wymiarowych
R
N
:
kxk
p
=
Ã
N
X
n=1
|x
n
|
p
!
1/p
,
dla
p = 1, 2, . . . , ∞
Wa˙zne
p
-normy
:
kxk
1
=
N
X
n=1
|x
n
|
- norma pierwsza,
kxk
2
=
v
u
u
t
N
X
n=1
|x
n
|
2
- norma Euklidesa,
kxk
∞
=
sup
n=1,2,...,N
|x
n
|
- norma maksimum (Czebyszewa)
Relacje pomi ˛edzy normami:
kxk
∞
6 kxk
2
6 kxk
1
6
√
N kxk
2
6 N kxk
∞
L.J. Opalski, “Wst ˛ep do metod numerycznych - materiały do wykładu”
wersja: 16.III.2010
Rozwi ˛
azywanie układów równa´n liniowych
24
Normy macierzy
A ∈ R
M ×N
indukowane
przez normy wektorów:
kAk
p
= sup
x6=0
kAxk
p
kxk
p
= max
kxk=1
kAxk
p
kAk
1
=
max
n=1,2,...,N
M
X
m=1
|a
mn
|
kAk
2
= max
x6=0
kAxk
2
kxk
2
=
r
max
i=1,...,N
λ
i
(A
t
A)
gdzie
λ
i
(B)
oznacza warto´s´c własn ˛
a macierzy
B
tzn. liczb ˛e spełniaj ˛
ac ˛
a równanie
Bv
i
= λv
i
dla pewnego wektora (własnego)
v
i
kAk
∞
=
max
m=1,2,...,M
N
X
n=1
|a
mn
|
Podstawowa relacja pomi ˛edzy normami:
kAk
2
2
6 kAk
1
· kAk
∞
Z definicji normy:
kAxk 6 kAk · kxk
O normie macierzy i normie wektora, spełniaj ˛
acych powy˙zszy warunek, mówimy ˙ze s ˛
a
zgodne
.
L.J. Opalski, “Wst ˛ep do metod numerycznych - materiały do wykładu”
wersja: 16.III.2010
Rozwi ˛
azywanie układów równa´n liniowych
25
Norma Frobeniusa (Euklidesa) macierzy:
kAk
E
=
v
u
u
t
M
X
m=1
N
X
n=1
|a
mn
|
2
jest zgodna
z norm ˛
a Euklidesa wektorów.
Własno´sci:
kAk
2
6 kAk
E
6
p
min(M, N ) · kAk
2
Uwaga:
Realizacja norm macierzy w programie
MATLAB
:
norma
wywołanie
sposób realizacji
A
1
norm(A,1)
max
(sum(abs(A)))
A
2
norm(A,2)
max
(svd(A))
A
E
norm(A,’fro’) sqrt
(sum(diag(A
0
∗ A)))
A
∞
norm(A,inf)
max
(sum(abs(A)))
L.J. Opalski, “Wst ˛ep do metod numerycznych - materiały do wykładu”
wersja: 16.III.2010
Rozwi ˛
azywanie układów równa´n liniowych
26
4.2.2 Niedokładno ´s ´c rozwi ˛
azania
•
Wpływ zaburze ´n wektora
b
na rozwi ˛
azanie
x
układu równa ´n
Ax = b
:
A(x + ∆x) = b + ∆b
⇒
A · ∆x = ∆b
⇒
∆x = A
−1
· ∆b
k∆xk 6 kA
−1
k · k∆bk,
oraz
kbk 6 kAk · kxk,
tzn.
1
kxk
6 kAk ·
1
kbk
St ˛
ad
k∆xk
kxk
6 cond(A)
k∆bk
kbk
cond(A) = kA
−1
k · kAk
-
wska´znik uwarunkowania
•
Dla
małych zaburze ´n
, tzn. gdy
kA
−1
k · k∆Ak < 1
, wpływ niedokładno´sci macierzy
A
mo˙zna
oszacowa´c nast ˛epuj ˛
aco:
k∆xk
kxk
6
cond(A)
k∆Ak
kAk
1 − cond(A)
k∆Ak
kAk
Uwaga: wska´znik uwarunkowania liczy w programie
MATLAB
funkcja
cond
L.J. Opalski, “Wst ˛ep do metod numerycznych - materiały do wykładu”
wersja: 16.III.2010
Rozwi ˛
azywanie układów równa´n liniowych
27
Przykład
Dany jest układ równa ´n
A · x = b
z macierz ˛
a Hilberta
A
, o elementach
a
i,j
= 1/(i + j − 1)
. Wektor
b
dobrano tak, by wektor rozwi ˛
aza ´n dokładnych (
x
) składał si ˛e z jedynek.
2
4
6
8
10
12
10
0
10
5
10
10
10
15
10
20
n
Propagacja niedokladnosci b
cond(A)
T(A)
Współczynnik propagacji niedokładno´sci wektora
b
szacowano wg:
T =
k˜
x − xk
kxk
/
k˜b − bk
kbk
wykorzystuj ˛
ac losowe zaburzenia wektora
b
(amplituda
wzgl ˛edna: 0.05%), ozn.
˜b
.
2
4
6
8
10
12
10
−20
10
−15
10
−10
10
−5
10
0
10
5
n
Maks. dokl. rozw.
||x1−x||
cond(A)*eps
Niedokładno´s´c rozwi ˛
azania numerycznego
x1
, a
zgrubne oszacowanie - za pomoc ˛
a:
cond(A) ∗ eps
.
L.J. Opalski, “Wst ˛ep do metod numerycznych - materiały do wykładu”
wersja: 16.III.2010
Rozwi ˛
azywanie układów równa´n liniowych
28
4.2.3 Poprawno ´s ´c numeryczna algorytmów sko ´nczonych
1. Algorytm Gaussa (oraz metoda LU) z przestawianiem równa ´n i/lub zmiennych s ˛
a numerycznie
poprawne, a równowa˙zne (bł ˛edom zaokr ˛
agle ´n operacji zmiennopozycyjnych) zaburzenie macierzy
układu spełnia warunek:
k∆Ak
1
kAk
1
6 O(n
3
) · eps
2. Algorytm rozkładu QR u˙zywaj ˛
acy przekształce ´n Hauseholdera jest numerycznie poprawny.
k∆Ak
E
kAk
E
6 cn
2
· eps,
gdzie
c ≈ 4
Dla ułatwienia dowodzenia poprawno´sci algorytmów przydatne jest poni˙zsze
Twierdzenie
Algorytm rozwi ˛
azywania układu
N
równa ´n liniowych jest numerycznie poprawny wtedy i tylko wtedy, gdy
istnieje stała
K
taka, ˙ze dla ka˙zdego układu równa ´n z tej klasy obliczony wektor reszt:
˜
r = fl(Ax − b)
spełnia warunek:
k˜
rk 6 K · kAk · k˜
xk · eps
, a
K = O(N
3
)
.
Zadanie (dla zainteresowanych). Pokaza´c, rozwi ˛
azuj ˛
ac układy równa ´n z macierz ˛
a Hilberta (z
poprzedniego przykładu), ˙ze algorytmy
MATLAB
a: LU, QR oraz algorytm obsługuj ˛
acy operator
\
zachowuj ˛
a si ˛e jak algorytmy poprawne numerycznie.
L.J. Opalski, “Wst ˛ep do metod numerycznych - materiały do wykładu”
wersja: 16.III.2010
Rozwi ˛
azywanie układów równa´n liniowych
29
4.3 Iteracyjne metody rozwi ˛
azywania układów równa ´n liniowych
4.3.1 Pochodzenie i własno ´sci wielkich zada ´n
Pochodzenie wielkich układów równa ´n liniowych
•
równania du˙zych sieci (np. elektrycznych)
•
dyskretyzacja równa ´n ró˙zniczkowych cz ˛
astkowych
Cechy wielkich równa ´n liniowych ułatwiaj ˛
ace rozwi ˛
azywanie:
•
małe wypełnienie (procent elementów niezerowych)
•
struktura blokowa lub pasmowa macierzy
•
dodatnia okre´slono´s´c:
∀x 6= 0 : x
T
Ax > 0
•
słaba diagonalna dominacja:
|a
mm
| >
P
N
n=1,n6=m
|a
mn
|
•
nieredukowalno´s´c macierzy. Macierz kwadratow ˛
a
A
nazywamy redukowaln ˛
a, je´sli istnieje macierz
permutacji
P
i macierze kwadratowe
B
,
C
takie, ˙ze
P
−1
AP =
·
B C
0 C
¸
•
znane jest dobre przybli˙zenie rozwi ˛
azania
L.J. Opalski, “Wst ˛ep do metod numerycznych - materiały do wykładu”
wersja: 16.III.2010
Rozwi ˛
azywanie układów równa´n liniowych
30
4.3.2 Wprowadzenie do algorytmów iteracyjnych
Przykład
Nale˙zy znale´z´c rozwi ˛
azanie
x
równania
a · x = 1
, gdzie
a
jest pewnym parametrem
Równanie to jest równowa˙zne nast ˛epuj ˛
acemu:
x = 1 + (1 − a) · x
Przyjmijmy jakie´s przybli˙zenie pocz ˛
atkowe rozwi ˛
azania, np.
x
(0)
= 0
;
kolejne przybli˙zenia wyznaczymy z rekurencyjnej formuły:
x
(n+1)
= 1 + (1 − a) · x
(n)
tzn.
x
(1)
= 1 + (1 − a) · x
(0)
= 1
,
x
(2)
= 1 + (1 − a) · x
(1)
= 2 − a
, ...
Iteracyjne poszukiwanie rozwi ˛
azania jest przydatne, je´sli
granica generowanego ci ˛
agu jest poszukiwanym rozwi ˛
azaniem
w niewielkiej liczbie iteracji mo˙zna uzyska´c dobre przybli˙zenie rozwi ˛
azania
.
L.J. Opalski, “Wst ˛ep do metod numerycznych - materiały do wykładu”
wersja: 16.III.2010
Rozwi ˛
azywanie układów równa´n liniowych
31
Zbie˙zno´s´c aperiodyczna
Zbie˙zno´s´c oscylacyjna
0
2
4
6
8
10
1
1.5
2
Numer iteracji (n)
x
(n)
dla rownania: 0.5 * x = 1
0
2
4
6
8
10
10
−4
10
−2
10
0
Numer iteracji (n)
abs(x
(n)
−2)
0
2
4
6
8
10
0.6
0.8
1
Numer iteracji (n)
x
(n)
dla rownania: 1.5 * x = 1
0
2
4
6
8
10
10
−4
10
−2
10
0
Numer iteracji (n)
abs(x
(n)
−0.666667)
Rozbie˙zno´s´c ci ˛
agu
0
2
4
6
8
10
−50
0
50
Numer iteracji (n)
x
(n)
dla rownania: 2.5 * x = 1
0
2
4
6
8
10
10
−2
10
0
10
2
Numer iteracji (n)
abs(x
(n)
−0.4)
Zachowanie ci ˛
agu
{x
(n)
}
dla kilku warto´sci parametru
a
L.J. Opalski, “Wst ˛ep do metod numerycznych - materiały do wykładu”
wersja: 16.III.2010
Rozwi ˛
azywanie układów równa´n liniowych
32
4.3.3 Ogólna idea iteracyjnego rozwi ˛
azywania układów równa ´n
Zadanie do rozwi ˛
azania:
Ax = b
Z macierzy
A
zostaje wydzielona pewna macierz
N
A = N + (A − N)
Nx
(i+1)
= b − (A − N)x
(i)
|
{z
}
d
(i)
x
(i+1)
= (I − N
−1
A)
|
{z
}
B
x
(i)
+ N
−1
b
| {z }
c
Warunki stosowalno´sci:
•
układ równa ´n
Nx
(i+1)
= d
(i)
jest łatwy (tani) do sformułowania i rozwi ˛
azania
•
ci ˛
ag przybli˙ze ´n
x
(i)
jest (i to dostatecznie szybko) zbie˙zny do rozwi ˛
azania:
lim
i→∞
x
(i)
= x
∞
= A
−1
b
L.J. Opalski, “Wst ˛ep do metod numerycznych - materiały do wykładu”
wersja: 16.III.2010
Rozwi ˛
azywanie układów równa´n liniowych
33
Twierdzenie
Liniowa metoda iteracyjna
x
(i+1)
= B · x
(i)
+ c
jest zbie˙zna do rozwi ˛
azania wtedy i tylko wtedy, gdy promie ´n spektralny:
ρ(B) < 1
oraz zachodzi warunek zgodno´sci:
ˆ
x = Bˆ
x + w
, gdzie
ˆ
x
jest rozwi ˛
azaniem dokładnym.
Uwagi:
• ρ(B) = {max(λ) | λ ∈ Spect(B)}
gdzie widmo
Spect
macierzy jest zbiorem jej
warto´sci własnych.
•
Szybko´s´c zbie˙zno´sci jest tym wi ˛eksza im promie ´n spektralny macierzy
B
jest mniejszy.
L.J. Opalski, “Wst ˛ep do metod numerycznych - materiały do wykładu”
wersja: 16.III.2010
Rozwi ˛
azywanie układów równa´n liniowych
34
4.3.4 Metoda Jacobiego
A = L + U + D
L
- macierz dolna trójk ˛
atna (z zerow ˛
a diagonal ˛
a)
U
- macierz górna trójk ˛
atna (z zerow ˛
a diagonal ˛
a)
D
- macierz diagonalna (nieosobliwa)
D · x
(i+1)
= −(L + U) · x
(i)
+ b
|
{z
}
d
(i)
albo
x
(i+1)
= −D
−1
(L + U)
|
{z
}
B
·x
(i)
+ D
−1
· b
| {z }
c
WKiW zbie˙zno´sci:
•
macierz
A
jest nieredukowalna i diagonalnie słabo dominuj ˛
aca, albo
•
macierz
A
jest silnie diagonalnie dominuj ˛
aca
L.J. Opalski, “Wst ˛ep do metod numerycznych - materiały do wykładu”
wersja: 16.III.2010
Rozwi ˛
azywanie układów równa´n liniowych
35
Przykład
Układ równa ´n do rozwi ˛
azania:
3
−2
0
0
−2
5
−2
0
0 −2
5
−2
0
0 −2
2
|
{z
}
A=L+D+U
x
1
x
2
x
3
x
4
| {z }
x
=
0
0
0
1
| {z }
b
Równania iteracji metody Jacobiego:
3
0 0 0
0
5
0 0
0 0
5
0
0 0 0
2
|
{z
}
D
x
(i+1)
= −
0 −2
0
0
−2
0 −2
0
0 −2
0 −2
0
0 −2
0
|
{z
}
L+U
x
(i)
+
0
0
0
1
| {z }
b
L.J. Opalski, “Wst ˛ep do metod numerycznych - materiały do wykładu”
wersja: 16.III.2010
Rozwi ˛
azywanie układów równa´n liniowych
36
4.3.5 Metoda Gaussa-Seidla
(L + D) · x
(i+1)
= −U · x
(i)
+ b
|
{z
}
d
(i)
albo
x
(i+1)
= −(L + D)
(−1)
· U
|
{z
}
B
·x
(i)
+ (L + D)
(−1)
· b
|
{z
}
c
WKiW zbie˙zno´sci:
•
macierz
A
jest dodatnio okre´slona, albo
•
macierz
A
jest silnie diagonalnie dominuj ˛
aca
L.J. Opalski, “Wst ˛ep do metod numerycznych - materiały do wykładu”
wersja: 16.III.2010
Rozwi ˛
azywanie układów równa´n liniowych
37
Przykład
Układ równa ´n do rozwi ˛
azania:
3
−2
0
0
−2
5
−2
0
0 −2
5
−2
0
0 −2
2
|
{z
}
A=L+D+U
x
1
x
2
x
3
x
4
| {z }
x
=
0
0
0
1
| {z }
b
Równania iteracji metody Gaussa-Seidla:
3
0
0 0
−2
5
0 0
0 −2
5
0
0
0 −2 2
|
{z
}
L+D
x
(i+1)
= −
0 −2
0
0
0
0 −2
0
0
0
0 −2
0
0
0
0
|
{z
}
U
x
(i)
+
0
0
0
1
| {z }
b
L.J. Opalski, “Wst ˛ep do metod numerycznych - materiały do wykładu”
wersja: 16.III.2010
Rozwi ˛
azywanie układów równa´n liniowych
38
4.3.6 Metody relaksacji
N(ω) =
1
ω
(D + ωL)
gdzie parametr relaksacji
ω ∈ (0, 2)
.
Równania iteracji metod relaksacji:
Nx
(i+1)
= b − (A − N)x
(i)
= b − (U +
ω − 1
ω
· D)x
(i)
albo
x
(i+1)
= (I − N
−1
A)
|
{z
}
B
x
(i)
+ N
−1
b
| {z }
c
Metoda kolejnych nadrelaksacji:
ω > 1
L.J. Opalski, “Wst ˛ep do metod numerycznych - materiały do wykładu”
wersja: 16.III.2010
Rozwi ˛
azywanie układów równa´n liniowych
39
Przykład
Układ równa ´n do rozwi ˛
azania:
3 −2
0
0
−2
5 −2
0
0 −2
5 −2
0
0 −2
2
|
{z
}
A=L+D+U
x
1
x
2
x
3
x
4
| {z }
x
=
0
0
0
1
| {z }
b
Równania iteracji metod relaksacji:
3/ω
0
0
0
−2 5/ω
0
0
0
−2 5/ω
0
0
0
−2 2/ω
|
{z
}
N=L+
1
ω
D
x
(i+1)
= −
3
ω−1
ω
−2
0
0
0 5
ω−1
ω
−2
0
0
0 5
ω−1
ω
−2
0
0
0 2
ω−1
ω
|
{z
}
U+
ω−1
ω
D
x
(i)
+
0
0
0
1
| {z }
b
L.J. Opalski, “Wst ˛ep do metod numerycznych - materiały do wykładu”
wersja: 16.III.2010
Rozwi ˛
azywanie układów równa´n liniowych
40
Porównanie szybko ´sci zbie˙zno ´sci algorytmów iteracyjnych
0
20
40
60
80
100
10
−20
10
−15
10
−10
10
−5
10
0
Liczba iteracji n
||x
(n)
−x||
Jacobi
Gauss−Seidel
SOR,
ω
=1.2
L.J. Opalski, “Wst ˛ep do metod numerycznych - materiały do wykładu”
wersja: 16.III.2010
Document Outline
- Rozwiazywanie układów równan liniowych
Wyszukiwarka
Podobne podstrony:
LAB 2 zad domowe, Lab2
GIiZK 0809 przydzial tematow zad domowego
METROL Zad Domowe2009 10 gr28
zad domowe-polecenie-wyklad, Magiczny Plik, 5 semestr, Rynki finansowe, zad dom
zad domowe
METROL Zad Domowe2008-09Tabele, simr pw, II rok, Metrologia, Prace metrologia, Pomoce projekty metro
Fizyka+lab+zad+12 (1)
METROL Zad Domowe2009 10 gr23
METROL Zad Domowe2009-10 gr28, simr pw, II rok, Metrologia, Prace metrologia, Pomoce projekty metrol
Ćwiczenia z MS Word Lab zad 3, Technologia Informacyjna semestr 1 oraz Informatyka i komputerowe ws
METROL Zad Domowe2009-10 gr23, SIMR, 3 semestr, metrologia, metrologia
zad domowe nr 1
METROL Zad Domowe2009-10 gr28
1 Zad domowe dodatkowe
Zad domowe III I IV I& LISTOPADA
Zad domowe 9, Studia, Rok I, Informatyka, semestr I
materiałoznawstwo zad domowe
rachunek, Lab 1 zad
Zad domowe korelacje
więcej podobnych podstron