Instytut Automatyki Politechniki Łódzkiej - Metody Numeryczne w Inżynierii wykład 2
W2 - 1
Rozwiązywanie układów równań liniowych
⎪
⎪
⎩
⎪⎪
⎨
⎧
=
+
+
+
=
+
+
+
=
+
+
+
n
n
nn
n
n
n
n
n
n
b
x
a
x
a
x
a
b
x
a
x
a
x
a
b
x
a
x
a
x
a
"
"
"
"
2
2
1
1
2
2
2
22
1
21
1
1
2
12
1
11
0
2
1
2
1
2
1
2
22
21
1
12
11
≠
=
⎥
⎥
⎥
⎥
⎦
⎤
⎢
⎢
⎢
⎢
⎣
⎡
=
⎥
⎥
⎥
⎥
⎦
⎤
⎢
⎢
⎢
⎢
⎣
⎡
⎥
⎥
⎥
⎥
⎦
⎤
⎢
⎢
⎢
⎢
⎣
⎡
A
det
,
b
Ax
b
b
b
x
x
x
a
a
a
a
a
a
a
a
a
n
n
nn
n
n
n
n
"
"
"
"
"
"
"
"
"
Instytut Automatyki Politechniki Łódzkiej - Metody Numeryczne w Inżynierii wykład 2
W2 - 2
Układy trójkątne
⎥
⎥
⎥
⎥
⎦
⎤
⎢
⎢
⎢
⎢
⎣
⎡
=
⎥
⎥
⎥
⎥
⎦
⎤
⎢
⎢
⎢
⎢
⎣
⎡
⎥
⎥
⎥
⎥
⎦
⎤
⎢
⎢
⎢
⎢
⎣
⎡
n
n
nn
n
n
b
b
b
x
x
x
u
u
u
u
u
u
"
"
"
"
"
"
"
"
"
2
1
2
1
2
22
1
12
11
0
0
0
1
1
1
1
,...,
n
,
n
i
,
x
u
b
u
x
n
i
k
k
ik
i
ii
i
−
=
⎥⎦
⎤
⎢⎣
⎡ −
=
∑
+
=
⎥
⎥
⎥
⎥
⎦
⎤
⎢
⎢
⎢
⎢
⎣
⎡
=
⎥
⎥
⎥
⎥
⎦
⎤
⎢
⎢
⎢
⎢
⎣
⎡
⎥
⎥
⎥
⎥
⎦
⎤
⎢
⎢
⎢
⎢
⎣
⎡
n
n
nn
n
n
b
b
b
x
x
x
l
l
l
l
l
l
"
"
"
"
"
"
"
"
"
2
1
2
1
2
1
22
21
11
0
0
0
n
,...,
,
i
,
x
l
b
l
x
i
k
k
ik
i
ii
i
2
1
1
1
1
=
⎥⎦
⎤
⎢⎣
⎡ −
=
∑
−
=
Instytut Automatyki Politechniki Łódzkiej - Metody Numeryczne w Inżynierii wykład 2
W2 - 3
Schemat eliminacji Gaussa z częściowym wyborem elementu
głównego
[
]
b
A
A
)
(
=
1
Postępowanie w k-tym kroku
)
k
(
)
k
(
A
A
1
+
→
1. Wybrać r:
)
k
(
ik
n
i
k
)
k
(
rk
a
max
a
≤
≤
=
2. Przestawić wiersze k i r , przestawienie zapamiętać
3. Obliczyć
n
,...,
k
,
k
i
,
a
a
m
)
k
(
kk
)
k
(
ik
k
,
i
2
1
+
+
=
=
4. Obliczyć
1
2
1
2
1
1
+
+
+
=
+
+
=
−
=
+
n
,...,
k
,
k
j
n
,...,
k
,
k
i
,
a
m
a
a
)
k
(
kj
ik
)
k
(
ij
)
k
(
ij
Instytut Automatyki Politechniki Łódzkiej - Metody Numeryczne w Inżynierii wykład 2
W2 - 4
Ostatecznie
⎥
⎥
⎥
⎥
⎦
⎤
⎢
⎢
⎢
⎢
⎣
⎡
=
⎥
⎥
⎥
⎥
⎦
⎤
⎢
⎢
⎢
⎢
⎣
⎡
⎥
⎥
⎥
⎥
⎦
⎤
⎢
⎢
⎢
⎢
⎣
⎡
)
n
(
n
)
(
)
(
n
)
n
(
nn
)
(
n
)
(
)
(
n
)
(
)
(
b
b
b
x
x
x
a
a
a
a
a
a
"
"
"
"
"
"
"
"
"
2
2
1
1
2
1
2
2
2
21
1
1
1
12
1
11
0
0
0
⎥
⎥
⎥
⎥
⎥
⎥
⎦
⎤
⎢
⎢
⎢
⎢
⎢
⎢
⎣
⎡
=
−
1
0
1
0
0
1
0
0
0
1
1
2
1
32
31
21
nn
n
n
m
m
m
m
m
m
L
"
#
%
#
#
#
"
"
"
Instytut Automatyki Politechniki Łódzkiej - Metody Numeryczne w Inżynierii wykład 2
W2 - 5
Sumy kontrolne:
∑
+
=
=
1
1
)
1
(
)
1
(
:
n
j
ij
i
a
s
,
n
i
,...,
3
,
2
,
1
=
[
]
s
b
A
A
=
)
1
(
2
,...,
3
,
2
,
,...,
3
,
2
,
)
1
(
1
1
)
1
(
)
2
(
+
=
=
−
=
n
j
n
i
a
m
a
a
j
i
ij
ij
n
i
s
m
s
s
i
i
i
,...,
3
,
2
,
)
1
(
1
1
)
1
(
)
2
(
=
−
=
[
]
n
i
a
a
m
a
a
m
a
s
n
j
ij
n
j
j
i
ij
n
j
j
i
n
j
ij
i
,...,
3
,
2
,
1
1
)
2
(
1
1
)
1
(
1
1
)
1
(
1
1
)
1
(
1
1
1
1
)
1
(
)
2
(
=
=
−
=
−
=
∑
∑
∑
∑
+
=
+
=
+
=
+
=
Instytut Automatyki Politechniki Łódzkiej - Metody Numeryczne w Inżynierii wykład 2
W2 - 6
Rozkład trójkątny
⎥
⎥
⎥
⎥
⎦
⎤
⎢
⎢
⎢
⎢
⎣
⎡
=
⎥
⎥
⎥
⎥
⎦
⎤
⎢
⎢
⎢
⎢
⎣
⎡
=
kk
k
k
k
k
k
nn
n
n
n
n
a
a
a
a
a
a
a
a
a
A
,
a
a
a
a
a
a
a
a
a
A
"
"
"
"
"
"
"
"
"
"
"
"
"
"
2
1
2
22
21
1
12
11
2
1
2
22
21
1
12
11
Twierdzenie:
Jeśli 1
2
1
0
−
=
≠
n
,...,
,
k
,
A
det
k
, to istnieje dokładnie jeden rozkład
LU
A
=
taki, że L jest macierzą trójkątna dolną z jedynkami na przekątnej ,
a U jest macierzą trójkątna górną.
Instytut Automatyki Politechniki Łódzkiej - Metody Numeryczne w Inżynierii wykład 2
W2 - 7
Macierzowy zapis eliminacji Gaussa bez wyboru elementu głównego
A
A
)
(
=
1
,
)
i
(
i
)
i
(
A
L
A
=
+1
i
i
m
m
L
i
,
n
i
,
i
i
#
"
"
"
#
%
#
#
#
"
"
"
"
#
#
#
%
#
"
"
⎥
⎥
⎥
⎥
⎥
⎥
⎥
⎦
⎤
⎢
⎢
⎢
⎢
⎢
⎢
⎢
⎣
⎡
−
−
=
+
1
0
0
0
1
0
0
0
1
0
0
0
0
1
1
i
i
m
m
L
i
n
i
i
i
#
"
"
"
#
%
#
#
#
"
"
"
"
#
#
#
%
#
"
"
⎥
⎥
⎥
⎥
⎥
⎥
⎥
⎥
⎦
⎤
⎢
⎢
⎢
⎢
⎢
⎢
⎢
⎢
⎣
⎡
=
+
−
1
0
0
0
1
0
0
0
1
0
0
0
0
1
,
,
1
1
LU
U
L
~
A
A
L
~
A
L
L
L
U
n
n
=
=
⇒
=
=
−
−
−
1
1
2
1
"
Jeżeli w trakcie eliminacji Gaussa wykonano przestawienia wierszy,
to
A
~
LU
=
gdzie
A
~
powstała z A przez przestawienia tych samych wierszy.
Instytut Automatyki Politechniki Łódzkiej - Metody Numeryczne w Inżynierii wykład 2
W2 - 8
Zastosowania rozkładu trójkątnego:
Rozwiązywanie układu równań liniowych Ax=b
[
]
b
A
A
)
(
=
1
→ El. Gaussa → 1trójkątny układ r-nań → x
Inaczej:
1) A
→ El. Gaussa, wyznaczenie macierzy L i U zapamiętanie
przestawień wierszy
2) wykonanie takich samych przestawień wierszy wektora b
b
b
~
→
3)
b
LUx
~
= , czyli
b
Ly
~
=
i
y
Ux
= 2 trójkątne układy równań do
rozwiązania
I
2
3
2
2
3
5
.
1
3
1
2
1
3
1
n
n
n
n
n
+
=
+
⎟
⎠
⎞
⎜
⎝
⎛
+
II
2
3
2
3
3
1
2
1
2
3
1
n
n
n
n
+
=
+
⎟
⎠
⎞
⎜
⎝
⎛
Instytut Automatyki Politechniki Łódzkiej - Metody Numeryczne w Inżynierii wykład 2
W2 - 9
Obliczanie wyznacznika
Jeżeli s to liczba wykonanych przestawień wierszy to:
)
(
)
2
(
22
)
1
(
11
)
1
(
)
det(
)
1
(
~
det
)
1
(
det
n
nn
s
s
s
a
a
a
LU
A
A
"
−
=
−
=
−
=
Instytut Automatyki Politechniki Łódzkiej - Metody Numeryczne w Inżynierii wykład 2
W2 - 10
Odwracanie macierzy:
Sposób 1:
[
]
1
2
1
:
−
=
=
A
x
x
x
X
n
#
"
#
#
[
]
⎥
⎥
⎥
⎥
⎦
⎤
⎢
⎢
⎢
⎢
⎣
⎡
=
=
1
0
0
0
1
0
0
0
1
:
2
1
"
%
"
"
#
"
#
#
n
e
e
e
I
j
j
e
Ax
=
[
]
I
A
→ el. Gaussa → n trójkatnych układów r-nań →
1
−
A
2
2
3
2
1
)
1
(
2
1
3
1
n
n
n
n
n
+
+
+
Instytut Automatyki Politechniki Łódzkiej - Metody Numeryczne w Inżynierii wykład 2
W2 - 11
Sposób 2:
1
1
1
1
)
(
−
−
−
−
=
=
L
U
LU
A
n
j
j
i
L
m
m
L
i
j
k
kj
ik
ij
ii
ij
,...,
1
,
)
(
1
)
(
1
1
1
+
=
⎥
⎦
⎤
⎢
⎣
⎡
−
=
∑
−
=
−
−
δ
1
,...,
1
,
)
(
)
(
1
)
(
1
1
1
−
=
⎥
⎦
⎤
⎢
⎣
⎡
−
=
∑
−
=
−
−
j
j
i
U
U
u
U
i
j
k
kj
ik
ij
ii
ij
δ
⎩
⎨
⎧
≠
=
=
j
i
j
i
j
i
0
1
,
δ
LU
A
=
~
→
1
1
1
1
)
(
~
−
−
−
−
=
=
L
U
LU
A
1
−
A
powstaje z
1
~
−
A
przez analogiczne do wykonanego przestawienia
wierszy przestawienie kolumn
2
3
3
1
n
n
+
Instytut Automatyki Politechniki Łódzkiej - Metody Numeryczne w Inżynierii wykład 2
W2 - 12
Błędy rozwiązań układu równań liniowych
Norma macierzy indukowana przez normę wektorową:
x
Ax
sup
A
x 0
≠
=
np.
j
,
i
j
,
i
i
i
a
max
A
,
x
max
x
=
=
Niech
b
Ax
= ,
b
)
x
x
)(
A
A
(
=
+
+
δ
δ
. Wtedy:
)
x
x
(
A
A
x
δ
δ
δ
+
−
=
−1
)
x
x
(
A
A
x
δ
δ
δ
+
≤
−1
A
A
A
A
)
x
x
(
x
δ
δ
δ
1
−
≤
+
Niech
b
Ax
= ,
b
b
)
x
x
(
A
δ
δ
+
=
+
. Wtedy:
b
x
A
δ
δ
=
b
A
x
δ
δ
1
−
≤
i
x
A
b
≤
b
b
A
A
x
x
δ
δ
1
−
≤
Wskaźnik uwarunkowania:
1
−
=
A
A
:
)
A
(
cond
Instytut Automatyki Politechniki Łódzkiej - Metody Numeryczne w Inżynierii wykład 2
W2 - 13
1
0.
)
A
(
cond
n
≤
ε
- dobre uwarunkowanie.
A = 1.2969 0.8648
0.2161 0.1441
b = 0.8642
0.1440
x=A\b
x = 2.0000
-2.0000
r=b-A*x
r =
1.0e-015 *
0.1110
0.0278
b1=b-[1e-8; -1e-8]
x1=A\b1
x1 =
0.9911
-0.4870
r1=b-A*x1
r1 =
1.0e-008 *
1.0000
-1.0000
b2=b+[1e-8; -1e-8]
x2=A\b2
x2 =
3.0089
-3.5130
r2=b-A*x2
r2 =
1.0e-007 *
-0.1000
0.1000
cond(A) = 2.4973e+008
inv(A) = 1.0e+008 * 0.1441 -0.8648
-0.2161 1.2969
Instytut Automatyki Politechniki Łódzkiej - Metody Numeryczne w Inżynierii wykład 2
W2 - 14
min(eig(A)) = 6.9396e-009, max(eig(A)) = 1.4410
2*1e-16*cond(A) = 4.9946e-008
Skalowanie
j
j
j
'
x
x
α
=
i i-te równanie mnożone przez
i
β
)
(
diag
D
),
(
diag
D
i
i
β
α
β
α
=
=
'
x
D
x
,
b
D
'
b
,
AD
D
'
A
,
'
b
'
x
'
A
α
β
α
β
=
=
=
=
Iteracyjne poprawianie rozwiązań:
)
i
(
x
i-te przybliżenie rozwiązania
Obliczyć z podwójną precyzją
)
i
(
)
i
(
Ax
b
r
−
=
Oszacować
1
1
10
10
−
−
≤
≤
−
>
>
)
i
(
)
i
(
k
)
i
(
j
n
j
k
r
max
Niech
)
i
(
k
v
)
i
(
10
będzie zaokrągleniem do pojedynczej precyzji
)
i
(
k
r
)
i
(
10
Obliczyć w pojedynczej precyzji
)
i
(
w
:
)
i
(
k
)
i
(
v
Aw
)
i
(
10
=
,
)
i
(
k
)
i
(
)
i
(
w
x
x
)
i
(
−
+
+
=
10
1
Wtedy
)
i
(
)
i
(
)
i
(
k
)
i
(
)
i
(
)
i
(
v
r
Aw
Ax
b
Ax
b
r
)
i
(
−
=
−
−
=
−
=
−
+
+
10
1
1
Instytut Automatyki Politechniki Łódzkiej - Metody Numeryczne w Inżynierii wykład 2
W2 - 15
Metody iteracyjne rozwiązywania układów r. liniowych-relaksacja:
1.
)
i
(
x
i-te przybliżenie rozwiązania
2. Obliczyć
)
i
(
)
i
(
Ax
b
r
−
=
3. Znaleźć składową o największym module:
)
i
(
l
r
, element o
największym module w l-tym wierszu macierzy A :
j
,
l
a
4. Obliczyć
)
i
(
x
1
+
zmieniając tylko j-tą składową w
)
i
(
x
:
j
,
l
)
i
(
l
)
i
(
j
)
i
(
j
a
r
x
x
+
=
+1
Wtedy
0
1
1
1
1
1
=
⇒
⎥
⎥
⎥
⎥
⎥
⎥
⎦
⎤
⎢
⎢
⎢
⎢
⎢
⎢
⎣
⎡
−
=
=
−
+
−
=
−
=
+
+
+
+
)
i
(
l
j
,
l
)
i
(
l
j
,
l
j
,
l
j
,
)
i
(
)
i
(
)
i
(
)
i
(
)
i
(
)
i
(
r
a
r
a
a
a
r
)
x
x
(
A
Ax
b
Ax
b
r
#
#