Instytut Automatyki Politechniki A
ódzkiej - Metody Numeryczne w Inżynierii wykład 2
Rozwiązywanie układów równań liniowych
a11x1 + a12 x2 + + a1n xn = b1
ż#
�#a x1 + a22 x2 + + a2n xn = b2
�#
21
�#

�#
�#an1x1 + an2 x2 + + ann xn = bn
�#
a11 a12 a1n x1 b1
Ą# ń#Ą# ń# Ą# ń#
ó#a
a22 a2n Ą#ó# x2 Ą# ó#b2 Ą#
21
ó# Ą#ó# Ą# ó# Ą#
=

ó# Ą#ó# Ą# ó# Ą#
ó#a
an2 ann Ą#ó#xn Ą# ó#bn Ą#
Ł# Ś#Ł# Ś# Ł# Ś#
n1
Ax = b, det A `" 0
W2 - 1
Instytut Automatyki Politechniki A
ódzkiej - Metody Numeryczne w Inżynierii wykład 2
Układy trójkątne
u11 u12 u1n x1 b1
Ą# ń#Ą# ń# Ą# ń#
ó#
0 u22 u2n Ą#ó# x2 Ą# ó#b2 Ą#
ó# Ą#ó# Ą# ó# Ą#
=

ó# Ą#ó# Ą# ó# Ą#
ó#
0 0 unn Ą#ó#xn Ą# ó#bn Ą#
Ł# Ś#Ł# Ś# Ł# Ś#
1
Ą#b - n
xi = i = n,n - 1,...,1
"u xk ń#,
uii ó# i k =i +1 ik Ą#
Ł# Ś#
l11 0 0 x1 b1
Ą# ń#Ą# ń# Ą# ń#
ó#l Ą#ó#
l22 0 x2 Ą# ó#b2 Ą#
21
ó# Ą#ó# Ą# ó# Ą#
=

ó# Ą#ó# Ą# ó# Ą#
ó#l
ln2 lnn Ą#ó#xn Ą# ó#bn Ą#
Ł# Ś#Ł# Ś# Ł# Ś#
n1
1
Ą#b - i-1 xk ń#,
xi = i = 1,2,...,n
"l
lii ó# i k =1 ik Ą#
Ł# Ś#
W2 - 2
Instytut Automatyki Politechniki A
ódzkiej - Metody Numeryczne w Inżynierii wykład 2
Schemat eliminacji Gaussa z częściowym wyborem elementu
głównego
A( 1 ) = [A b]
Postępowanie w k-tym kroku A( k ) A( k +1 )
1. Wybrać r:
( k (
ark ) = max aikk )
kd"id"n
2. Przestawić wiersze k i r , przestawienie zapamiętać
3. Obliczyć
(
aikk )
mi ,k = , i = k + 1,k + 2,...,n
( k
akk )
( ( (
aijk +1 ) = aijk ) - mikakjk ) ,
4. Obliczyć i = k + 1,k + 2,...,n
j = k + 1,k + 2,...,n + 1
W2 - 3
Instytut Automatyki Politechniki A
ódzkiej - Metody Numeryczne w Inżynierii wykład 2
Ostatecznie
( 1 ( 1 ( 1 (
x1
Ą# ń# Ą# ń#
a11 ) a12 ) a1n ) b1 1 )
Ą# ń#
ó#
( 2 ( 2 (
x2
0 a21 ) a2n ) Ą#ó# Ą# ó#b2 2 ) Ą#
ó# Ą# ó# Ą#
ó# Ą#
=
ó# Ą# ó# Ą#

ó# Ą#

ó# Ą# ó# Ą#
( n (
0 0 ann )Ś#ó#xn Ą# Ł#bnn )Ś#
Ł# Ś#
Ł#
1 0 0 0
Ą# ń#
ó#m
1 0 0Ą#
21
ó# Ą#
L = m31 m32 1 0
ó# Ą#
ó#
Ą#
ó# Ą#
ó# Ą#
mn2 mnn-1 1Ś#
Ł#m
n1
W2 - 4
Instytut Automatyki Politechniki A
ódzkiej - Metody Numeryczne w Inżynierii wykład 2
Sumy kontrolne:
n+1
(1) (1)
[]
si := i = 1,2,3,..., n A(1) = A b s
"a ,
ij
j=1
( ( (
aij2) = aij1) - mi1a11), i = 2,3,..., n, j = 2,3,..., n + 2
j
(2) (1) (
si = si - mi1s11), i = 2,3,..., n
n+1 n+1 n+1 n+1
(2) (1) (1) ( ( (2)
si = [ - mi1a11) =
aij1) ]
"a - mi1"a = " "a , i = 2,3,..., n
ij 1 j j ij
j=1 j=1 j=1 j=1
W2 - 5
Instytut Automatyki Politechniki A
ódzkiej - Metody Numeryczne w Inżynierii wykład 2
Rozkład trójkątny
a11 a12 a1n a11 a12 a1k
Ą# ń# Ą# ń#
ó#a ó#a
a22 a2n Ą# a22 a2k Ą#
21 21
ó# Ą#, Ak = ó# Ą#
A =

ó# Ą# ó# Ą#
ó#a ó#a
an2 ann Ą# ak 2 akk Ą#
Ł# Ś# Ł# Ś#
n1 k1
Twierdzenie:
Jeśli det Ak `" 0, k = 1,2,...,n - 1, to istnieje dokładnie jeden rozkład
A = LU
taki, że L jest macierzą trójkątna dolną z jedynkami na przekątnej ,
a U jest macierzą trójkątna górną.
W2 - 6
Instytut Automatyki Politechniki A
ódzkiej - Metody Numeryczne w Inżynierii wykład 2
Macierzowy zapis eliminacji Gaussa bez wyboru elementu głównego
A( 1 ) = A, A( i +1 ) = Li A( i )
1 0 0 0
Ą# ń#
ó#
Ą#
ó# Ą#
0 1 0 0
ó# Ą# i
1 0 0 0
Ą# ń#
ó#0 - mi+1,i 1 0Ą#
ó#
Ą#
ó# Ą#
Li =
ó# Ą#
ó#
Ą#
ó# Ą#
0 1 0 0
i
ó#0 - mn 0 1Ą#
ó#0 mi
1 0Ą#
,i
Ł# Ś#
+1,i
ó# Ą#
L-1 =
i

ó#
Ą#
ó# Ą#
i
ó#0 mn,i 0 1Ą#
Ł# Ś#

i
~ ~
U = Ln-1Ln-2 L1 A = LA �! A = L-1U = LU
Jeżeli w trakcie eliminacji Gaussa wykonano przestawienia wierszy,
to
~
LU = A
~
gdzie A powstała z A przez przestawienia tych samych wierszy.
W2 - 7
Instytut Automatyki Politechniki A
ódzkiej - Metody Numeryczne w Inżynierii wykład 2
Zastosowania rozkładu trójkątnego:
Rozwiązywanie układu równań liniowych Ax=b
A( 1 ) = [A b]
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
~
~, czyli Ly = b i Ux = y 2 trójkątne układy równań do
3) LUx = b
rozwiązania
1 1 1
�#
1 1 1
�#
n3 ś# + 2 n2 = n3 + n2
ś# ź#
n3 + n2 ś# + n2 = n3 +1.5n2 II
ś# ź#
I
3 2 3
3 2 3 �# #
�# #
W2 - 8
Instytut Automatyki Politechniki A
ódzkiej - Metody Numeryczne w Inżynierii wykład 2
Obliczanie wyznacznika
Jeżeli s to liczba wykonanych przestawień wierszy to:
~
(1 (2 (n
det A = (-1)s det A = (-1)s det(LU ) = (-1)s a11)a22) ann)
W2 - 9
Instytut Automatyki Politechniki A
ódzkiej - Metody Numeryczne w Inżynierii wykład 2
Odwracanie macierzy:
Sposób 1:
-1
[ ]
X = x1 x x := A
2 n
1 0 0
Ą# ń#
ó#0 1 0Ą#
ó# Ą#
[]
I = e1 e e :=
2 n

ó# Ą#
ó#0 0 1 Ą#
Ł# Ś#
Ax = ej
j
[ ]
A I el. Gaussa n trójkatnych układów r-nań A-1
1 1 1
n3 + n2(n +1) + n n2
3 2 2
W2 - 10
Instytut Automatyki Politechniki A
ódzkiej - Metody Numeryczne w Inżynierii wykład 2
Sposób 2:
-1
A-1 = (LU )-1 = U L-1
i-1
Ą# ń#
1
(L-1)ij = -
"m (L-1)kj Ś# i = j, j + 1,..., n
Ą#
mii ó#�ij k= j ik
Ł#
i-1
Ą# ń#
1
-1 -1
(U )ij = -
"(U )ik (U )kj Ś# i = j, j - 1,...,1
Ą#
uii ó#�ij k= j
Ł#
1 i = j
ż#
�i , j =
�#
�#0 i `" j
~ ~
-1
A = LU A-1 = (LU )-1 = U L-1
~
A-1 powstaje z A-1 przez analogiczne do wykonanego przestawienia
wierszy przestawienie kolumn
1
n3 + n2
3
W2 - 11
Instytut Automatyki Politechniki A
ódzkiej - Metody Numeryczne w Inżynierii wykład 2
Błędy rozwiązań układu równań liniowych
Norma macierzy indukowana przez normę wektorową:
Ax
A = sup
x
x `"0
np. x = max xi , A = max ai , j
i i , j
Niech Ax = b, ( A + �A )( x + �x ) = b. Wtedy:
�x = - A-1�A( x + �x ) �x d" A-1 �A ( x + �x )
�x �A
d" A A-1
( x + �x ) A
Niech Ax = b, A( x + �x ) = b + �b. Wtedy: A�x = �b
�x d" A-1 �b i b d" A x
�x �b
d" A A-1
x b
Wskaz
nik uwarunkowania: cond( A ) := A A-1
W2 - 12
Instytut Automatyki Politechniki A
ódzkiej - Metody Numeryczne w Inżynierii wykład 2
n� cond( A ) d" 0.1 - dobre uwarunkowanie.
A = 1.2969 0.8648 b1=b-[1e-8; -1e-8] b2=b+[1e-8; -1e-8]
0.2161 0.1441
b = 0.8642
0.1440 x1=A\b1 x2=A\b2
x=A\b x1 = x2 =
x = 2.0000 0.9911 3.0089
-2.0000 -0.4870 -3.5130
r=b-A*x r1=b-A*x1 r2=b-A*x2
r = r1 = r2 =
1.0e-015 * 1.0e-008 * 1.0e-007 *
0.1110 1.0000 -0.1000
0.0278 -1.0000 0.1000
cond(A) = 2.4973e+008
inv(A) = 1.0e+008 * 0.1441 -0.8648
-0.2161 1.2969
W2 - 13
Instytut Automatyki Politechniki A
ódzkiej - Metody Numeryczne w Inżynierii wykład 2
min(eig(A)) = 6.9396e-009, max(eig(A)) = 1.4410
2*1e-16*cond(A) = 4.9946e-008
Skalowanie
x = ą x' i i-te równanie mnożone przez �i
j j j
Dą = diag(ąi ), D� = diag( �i )
A' x' = b' , A' = D� ADą , b' = D� b, x = Dą x'
Iteracyjne poprawianie rozwiązań:
x( i ) i-te przybliżenie rozwiązania
Obliczyć z podwójną precyzją r( i ) = b - Ax( i )
( i ) ( i )
Oszacować 10-k > max rj( i ) > 10-k -1
1d" jd"n
( i )
Niech 10k v( i ) będzie zaokrągleniem do pojedynczej precyzji
( i )
10k r( i )
( i )
Obliczyć w pojedynczej precyzji w( i ): Aw( i ) = 10k v( i ),
( i )
x( i +1 ) = x( i ) + 10-k w( i )
( i )
( i )
Wtedy r( i +1 ) = b - Ax( i +1 ) = b - Ax( i ) - 10-k Aw( i ) = r - v( i )
W2 - 14
Instytut Automatyki Politechniki A
ódzkiej - Metody Numeryczne w Inżynierii wykład 2
Metody iteracyjne rozwiązywania układów r. liniowych-relaksacja:
1. x( i ) i-te przybliżenie rozwiązania
2. Obliczyć r( i ) = b - Ax( i )
3. Znalez
ć składową o największym module: rl( i ), element o
największym module w l-tym wierszu macierzy A : al , j
4. Obliczyć x( i +1 ) zmieniając tylko j-tą składową w x( i ):
rl( i )
x(j i +1 ) = x(j i ) +
al , j
Wtedy
r( i +1 ) = b - Ax( i +1 ) = b - Ax( i ) + A( x( i ) - x( i +1 ) ) =
a1, j
Ą# ń#
ó# Ą#

i )
ó# Ą#
al , j rl(
= r( i ) - ó# Ą# �! rl( i +1 ) = 0
al , j
ó# Ą#

ó# Ą#
ó# Ą#
l , j
Ł#a Ś#
W2 - 15