LAB 2 zad domowe WNUM W04

Rozwiązywanie układów równań liniowych 1
4 Rozwiązywanie układów równań liniowych
Standardowa postać układu równań:
�ł łł �ł łł �ł łł
A1,1 A1,2 � � � A1,n-1 An,n x1 b1
�ł śł �ł śł �ł śł
�ł
A2,1 A2,2 � � � A2,n-1 An,n śł �ł x2 śł �ł b2 śł
�ł śł �ł śł �ł śł
�ł śł �ł śł �ł śł
. . . . . . .
�ł śł �ł śł �ł śł
. . . . . . .
� =
. . . . . . .
�ł śł �ł śł �ł śł
�ł śł �ł śł �ł śł
�ł
An-1,1 An-1,2 � � � An-1,n-1 An-,n śł �ł xn-1 śł �ł bm-1 śł
�ł �ł �ł �ł �ł �ł
An,1 An,2 � � � An,n-1 An,n xn bm
A x
b
Z algebry wiadomo, że jeśli macierz A jest nieosobliwa, to jedyne rozwiązanie tego układu
wyraża się jako: x = A-1 � b , gdzie A-1 jest macierzą odwrotną do macierzy A.
L.J. Opalski, Wstęp do metod numerycznych - materiały do wykładu wersja: 16.III.2010
Rozwiązywanie układów równań liniowych 2
Wyznaczanie rozwiązania układu równań x = A-1 � b w dwóch krokach:
1. oblicz macierz C odwrotną do A, t.j.: C = A-1
2. oblicz rozwiązanie: x = C � b
nie jest korzystne ze względu na nakłady obliczeniowe i dokładność.
Uzasadnienie-przykład
Równanie 7 � x = 21 ma rozwiązanie całkowite: x = 21/7 = 3. Wyznaczanie rozwiązania
przy użyciu odwrotności współczynnika przy x, tzn.:
x = 7-1 � 21 H" 0.142857 � 21 = 2.99997, jest mniej dokładne i wymaga więcej operacji
arytmetycznych.
Dalej będzie mowa:
" o nieiteracyjnych algorytmach efektywnego i dokładnego rozwiązywania równań liniowych
" o tym jak na dokładność wyznaczania rozwiązania wpływa dokładność danych (elementów
macierzy A i wektora b) oraz wybór algorytmu.
" o algorytmach iteracyjnych i obszarach ich zastosowań
L.J. Opalski, Wstęp do metod numerycznych - materiały do wykładu wersja: 16.III.2010
Rozwiązywanie układów równań liniowych 3
4.1 Nieiteracyjne rozwiązywanie układów równań
4.1.1 Metoda eliminacji Gaussa
Eliminacja Gaussa jest algorytmem dwufazowym:
" Pierwsza faza przekształca obydwie strony układu równań tak, by uzyskać równoważny
układ z macierzą górną trójkątną (U):
U � x = z
W każdym kroku przekształceń zerowane są elementy kolejnej kolumny macierzy układu
poniżej diagonali
" W fazie drugiej rozwiązywany jest układ równań U � x = z względem kolejnych
zmiennych - począwszy od ostatniej.
L.J. Opalski, Wstęp do metod numerycznych - materiały do wykładu wersja: 16.III.2010
Rozwiązywanie układów równań liniowych 4
Przykład
�ł łł �ł łł �ł łł �ł łł
R1 : 3 1 6 x1 2 R1 bez zmian 1 0 0
�ł śł �ł śł �ł śł �ł śł
�ł śł �ł śł=�ł śł �ł śł
=�! =�!
R2 : 2 1 3 x2 7 R2 := R2 - R1 " 2 / 3 2/3 1 0
�ł �ł �ł �ł �ł �ł �ł �ł
R3 : 1 1 1 x3 4 R3 := R3 - R1 " 1 / 3 1/3 0 1
A x b
L(1)
�ł łł
R1 bez zmian
�ł łł �ł łł �ł łł
3 1 6
R1 : x1 2 1 0 0
�ł śł
R2 bez zmian
�ł śł �ł śł �ł śł
�ł śł
1
17 2
�ł śł=�ł śł �ł śł
0 -1
=�! =�!
�ł śł
R2 : x2 1 0
3
2
�ł �ł �ł �ł �ł �ł
3 3
�ł �ł
3
2 10 1
R3 := R3 - R2 "
R3 : 0 -1 x3 2 1
3 3 3
1
3
x
b(1) La"L(2)
A(1)
�ł łł
�ł łł �ł łł
3 1 6
R1 : x1 2 L(1)A(1) = A , L(1)b(1) = b
�ł śł
�ł śł �ł śł
�ł śł
1
17
�ł śł=�ł śł
Uwaga:
�ł 0 -1 śł
R2 : x2 L(2)A(2) = A , L(2)b(2) = b
3 �ł �ł �ł �ł
3
�ł �ł
R3 : x3 -8 LU = A , Lz = b
0 0 1
x
za"b(2)
Ua"A(2)
2- 1 �x2- 6 �x3
�ł łł
R1 : x1 =
3
19
�ł śł
17/3-( -1 )�x3
�ł śł
�! x =
-7
R2 : x2 =
�ł �ł
1/3
-8
-8
R3 : x3 =
1
L.J. Opalski, Wstęp do metod numerycznych - materiały do wykładu wersja: 16.III.2010
Rozwiązywanie układów równań liniowych 5
Przykład
Znajdzmy napięcia węzłowe U1, U2, U3, U4 w elektrycznej sieci liniowej:
Równania bilansu prądów (dla węzłow 1, 2, 3, 4)
dają:
G34 G23 2 G12 1
4 3
G1U1 + G12(U1 - U2) = 0
G3 G2 G1
I = 1
G2U2 + G12(U2 - U1) + G23(U2 - U3) = 0
G3U3 + G23(U3 - U2) + G3(U3 - U4) = 0
G4(U4 - U3) = 1
Przyjmijmy Gi = 1, Gij = 2, otrzymując
układ równań: Sieć elektryczną
�ł łł �ł łł �ł łł
3 -2 0 0 U1 0
4 2 3 2 2 2 1
�ł śł �ł śł �ł śł
�ł śł �ł śł �ł śł
-2 5 -2 0 U2 0
�ł śł �ł śł �ł śł
=
�ł śł �ł śł �ł śł
�ł śł �ł śł �ł śł
0 -2 5 -2 U3 0
�ł �ł �ł �ł �ł �ł
1 1 1 1
0 0 -2 2 U4 1
2�1 2
Geq,1 = =
2+1 3
L.J. Opalski, Wstęp do metod numerycznych - materiały do wykładu wersja: 16.III.2010
Rozwiązywanie układów równań liniowych 6
Etap I. Eliminacja zmiennych
1. Mnożąc równanie 1 przez 2/3 i dodając do drugiego:
�ł łł �ł łł �ł łł
4 2 3 2 2
3 -2 0 0 U1 0
�ł śł �ł śł �ł śł
�ł śł �ł śł �ł śł
0 11/3 -2 0 U2 0
�ł śł �ł śł �ł śł
= 5
�ł śł �ł śł �ł śł
1 1
�ł śł �ł śł �ł śł
3
0 -2 5 -2 U3 0
�ł �ł �ł �ł �ł �ł
0 0 -2 2 U4 1
2�5/3
10
Geq,2 = =
2+5/3 11
2. Mnożąc równanie 2 przez 6/11 i dodając do trzeciego:
�ł łł �ł łł �ł łł
4 2 3
3 -2 0 0 U1 0
�ł śł �ł śł �ł śł
�ł śł �ł śł �ł śł
0 11/3 -2 0 U2 0
�ł śł �ł śł �ł śł
= 21
�ł śł �ł śł �ł śł
1
�ł śł �ł śł �ł śł
11
0 0 43/11 -2 U3 0
�ł �ł �ł �ł �ł �ł
0 0 -2 2 U4 1
2�21/11
42
Geq,3 = =
2+21/11 43
L.J. Opalski, Wstęp do metod numerycznych - materiały do wykładu wersja: 16.III.2010
Rozwiązywanie układów równań liniowych 7
3. Mnożąc równanie 3 przez 22/43 i dodając do czwartego:
�ł łł �ł łł �ł łł
4
3 -2 0 0 U1 0
�ł śł �ł śł �ł śł
�ł śł �ł śł �ł śł
0 11/3 -2 0 U2 0
�ł śł �ł śł �ł śł
42
=
�ł śł �ł śł �ł śł
1
43
�ł śł �ł śł �ł śł
0 0 43/11 -2 U3 0
�ł �ł �ł �ł �ł �ł
0 0 0 42/43 U4 1
Etap II. Podstawienia odwrotne (wstecz)
�ł łł �ł łł �ł łł
3 -2 0 0 U1 0
�ł śł �ł śł �ł śł
�ł śł �ł śł �ł śł
0 11/3 -2 0 U2 0
�ł śł �ł śł �ł śł
=
�ł śł �ł śł �ł śł
�ł śł �ł śł �ł śł
0 0 43/11 -2 U3 0
�ł �ł �ł �ł �ł �ł
0 0 0 42/43 U4 1
Stąd:
42 43
U4 = 1 �! U4 =
43 42
43 11
U3 - 2U4 = 0 �! U3 =
11 21
11 6
U2 - 2U3 = 0 �! U2 =
3 21
4
3U1 - 2U2 = 0 �! U1 =
21
L.J. Opalski, Wstęp do metod numerycznych - materiały do wykładu wersja: 16.III.2010
Rozwiązywanie układów równań liniowych 8
Wykonalność algorytmu Gaussa
Przykład
�ł łł �ł łł �ł łł
1 2 3 6 1
�ł �ł �ł �ł �ł �ł
Układ równań 1 2 1 x = 4 ma dokładnie jedno rozwiązanie: 1
2 1 3 6 1
A b
Tymczasem po jednym kroku eliminacji Gaussa mamy:
�ł łł �ł łł
1 2 3 6
�ł �ł �ł �ł
0 0 -2 x = -2
0 -3 -3 -6
Drugiego kroku eliminacji nie można zrobić!
L.J. Opalski, Wstęp do metod numerycznych - materiały do wykładu wersja: 16.III.2010
Rozwiązywanie układów równań liniowych 9
Rozwiązanie problemu : przestawienie równań 2 i 3.
�ł łł �ł łł �ł łł �ł łł
1 0 0 1 2 3 1 0 0 6
�ł �ł �ł �ł �ł �ł �ł �ł
0 0 1 1 2 1 x = 0 0 1 4
0 1 0 2 1 3 0 1 0 6
P A P b
Am bm
gdzie P jest macierzą przestawień.
�ł łł �ł łł
1 2 3 6
�ł �ł �ł �ł
Am = 2 1 3 , bm = 6
1 2 1 4
Po jednym kroku eliminacji Gaussa mamy układ równań z macierzą górną trójkątną:
�ł łł �ł łł
1 2 3 6
�ł �ł �ł �ł
0 -3 -3 x = -6
0 0 -2 -2
U z
L.J. Opalski, Wstęp do metod numerycznych - materiały do wykładu wersja: 16.III.2010
Rozwiązywanie układów równań liniowych 10
Przestawienia równań wpływają nie tylko na wykonalność eliminacji Gaussa, ale i na dokładność
�ł łł �ł łł �ł łł
Przykład
10 -7 0 x1 7
�ł śł �ł śł �ł śł
�ł śł �ł śł=�ł śł
Rozwiązać układ równań korzystając z aryt-
-3 2.099 6 x2 3.901
�ł �ł �ł �ł �ł �ł
mometru wykonującego operacje zmiennopozy-
5 -1 5 x3 6
cyjne z 5 dziesiętnymi cyframi znaczącymi.
�ł łł �ł łł �ł łł
10 -7 0 x1 7
�ł śł �ł śł �ł śł
�ł śł �ł śł=�ł śł
Po pierwszym kroku eliminacji Gaussa:
0 -0.001 6 x2 6.001
�ł �ł �ł �ł �ł �ł
0 2.5 5 x3 2.5
�ł łł
�ł łł �ł łł
7
śł
10 -7 0 x1 �ł
�ł śł
�ł śł �ł śł
�ł śł �ł śł=�ł 6.001 śł
Po drugim kroku eliminacji Gaussa:
�ł śł
0 -0.001 6 x2 �ł
�ł �ł �ł �ł
śł
15006
�ł
0 0 15005 x3 �ł
z3
Podstawienia wstecz dają: x3 = 15006/15005 = 1.0001,
x2 = (6.001 - 6 " x3)/(-0.001) = -0.4, x1 = (7 - 7 � x2)/10 = 0.9800.
Wartości dokładne: x1 = 0, x2 = -1, x3 = 1.
Przyczyna problemu: zły wybór elementu centralnego w kroku 1 �! duża wrażliwość na błąd
zaokrągleń przy obliczaniu z3.
L.J. Opalski, Wstęp do metod numerycznych - materiały do wykładu wersja: 16.III.2010
Rozwiązywanie układów równań liniowych 11
�ł łł �ł łł �ł łł
10 -7 0 x1 7
�ł śł �ł śł �ł śł
�ł śł �ł śł=�ł śł
Postać układu równań po zamianie dru-
5 -1 5 x2 6
�ł �ł �ł �ł �ł �ł
giego i trzeciego równania
-3 2.099 6 x3 3.901
�ł łł �ł łł �ł łł
10 -7 0 x1 7
�ł śł �ł śł �ł śł
�ł śł �ł śł=�ł śł
Po pierwszym kroku eliminacji Gaussa:
0 2.5 5 x2 2.5
�ł �ł �ł �ł �ł �ł
0 -0.001 6 x3 6.001
�ł łł �ł łł �ł łł
10 -7 0 x1 7
�ł śł �ł śł �ł śł
�ł śł �ł śł=�ł śł
Po drugim kroku eliminacji Gaussa:
0 2.5 5 x2 2.5
�ł �ł �ł �ł �ł �ł
0 0 6.002 x3 6.002
Podstawienia wstecz dają wartości dokładne: x3 = 1, x2 = (2.5 - 5 � x3)/2.5 = -1,
x1 = (7 + 7 " x2)/10 = 0.
L.J. Opalski, Wstęp do metod numerycznych - materiały do wykładu wersja: 16.III.2010
Rozwiązywanie układów równań liniowych 12
Uwagi
" Przestawianie wierszy (bądz wierszy i kolumn), prowadzące do dużych (bezwzględnie)
wartości elementów centralnych (diagonalnych) zmniejsza przenoszenie błędów
zaokrągleń.
" W dwóch ważnych przypadkach można nie wykonywać przestawień:
jeżeli macierz główna A jest symetryczna i dodatnio określona, tzn.
AT = A, i xT Ax > 0, dla każdego x = 0

jeżeli macierz główna ma dominującą przekątną, tzn:
n
|ai,i| |ai,j|
j=1,j=i

L.J. Opalski, Wstęp do metod numerycznych - materiały do wykładu wersja: 16.III.2010
Rozwiązywanie układów równań liniowych 13
Nakłady obliczeniowe algorytmu eliminacji Gaussa
1 1
mnożenie i dzielenie: M = n3 + n2 - n
3 3
1 1 5
dodawanie i odejmowanie: D = n3 + n2 - n
3 2 6
Uwagi
" wzory Cramera wymagają (n + 1)! mnożeń
" MATLAB: rozwiązanie układu równań Ax = b za pomocą: x=A\b
L.J. Opalski, Wstęp do metod numerycznych - materiały do wykładu wersja: 16.III.2010
Rozwiązywanie układów równań liniowych 14
4.1.2 Wykorzystanie rozkładu LU
Załóżmy, że
A = L � U
gdzie:
L jest macierzą dolną trójkątną,
U jest macierzą górną trójkątną.
Rozwiązanie układu równań A � x = b uzyskuje się 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ęp do metod numerycznych - materiały do wykładu wersja: 16.III.2010
Rozwiązywanie układów równań liniowych 15
�ł łł �ł łł �ł łł �ł łł
1 1 1 x1 łł �ł 3 1 0 0 1 1 1
Przykład
�ł �ł �ł �ł �ł, A = �ł �ł �ł �ł
1 -1 1 x2 �ł = 1 1 1 0 � 0 -2 0
1 2 3 x3 6 1 -1/2 1 0 0 2
x
A b L U
�ł łł �ł łł �ł łł
1 0 0 1 1 1 x1 łł �ł 3
�ł �ł �ł �ł �ł �ł �ł,
1 1 0 � 0 -2 0 x2 �ł = 1
1 -1/2 1 0 0 2 x3 6
x
L U b
z
�ł łł �ł łł �ł łł �ł łł �ł łł
1 0 0 z1 3 z1 3
�ł �ł �ł �ł �ł �ł �ł �ł
1 1 0 � z2 �ł = 1 =�! z2 �ł = -2
1 -1/2 1 z3 6 z3 2
z z
L b
�ł łł �ł łł �ł łł �ł łł
1 1 1 x1 łł �ł 3 x1 1
�ł �ł �ł �ł �ł �ł �ł �ł
0 -2 0 x2 �ł = -2 =�! x2 �ł = 1
0 0 2 x3 2 x3 1
x z x
U
L.J. Opalski, Wstęp do metod numerycznych - materiały do wykładu wersja: 16.III.2010
Rozwiązywanie układów równań liniowych 16
Nakłady obliczeniowe dla metody LU
rozkład macierzy:
1 1
mnożenie i dzielenie: M = n3 - n
3 3
1 1 1
dodawanie i odejmowanie: D = n3 - n2 + n
3 2 6
podstawienia:
mnożenie i dzielenie: M = n2
dodawanie i odejmowanie: D = n2 - n
Uwaga
Nakłady obliczeniowe na rozwiązanie układu równań za pomocą metody LU są takie same jak
dla eliminacji Gaussa.
L.J. Opalski, Wstęp do metod numerycznych - materiały do wykładu wersja: 16.III.2010
Rozwiązywanie układów równań liniowych 17
Rozkład LU w programie MATLAB
1.[L,U,P]=lu(A)rozkłada macierz A tak, że PA = LU, gdzie L jest macierzą
trójkątną dolną, U jest macierzą trójkątną górną, a P jest macierzą przestawień.
Ax = b �! PA x = Pb �! L Ux = Pb
LU z
Wektor x wynika z rozwiązania dwóch układów z trójkątnymi macierzami:
Lz = Pb, Ux = z
2.[Lm,U]=lu(A)oblicza rozkład macierzy A na iloczyn A = LmU, przy czym U jest
macierzą trójkątną górną, a Lm jest iloczynem macierzy przestawień i macierzy dolnej
trójkątnej. Rozwiązanie układu równań Ax = b uzyskuje się z dwóch układów z
trójkątnymi macierzami:
Lmz = b, Ux = z
L.J. Opalski, Wstęp do metod numerycznych - materiały do wykładu wersja: 16.III.2010
Rozwiązywanie układów równań liniowych 18
Przykład�ł łł�ł łł
1 1 1 x1 łł �ł 3
A=[1 ,1 ,1 ; 1 , 2 ,3;1 , -1 ,1 ] ;
�ł �ł�ł �ł �ł
1 2 3 x2 �ł = 6
b = [ 3 ; 1 ; 6 ] ;
1 -1 1 x3 1
[ L , U, P]= lu (A)
[Lm,Um]= lu (A) % L =
% Lm = % 1.0000 0 0
% 1.0000 0 0 % 1.0000 1.0000 0
% 1.0000 -0.5000 1.0000 % 1.0000 -0.5000 1.0000
% 1.0000 1.0000 0 % U =
% Um = % 1 1 1
% 1 1 1 % 0 -2 0
% 0 -2 0 % 0 0 2
% 0 0 2 % P =
z2=Lm\ b ; % 1 0 0
x2=Um\ z2 % 0 0 1
% 1 % 0 1 0
% 1 z1=L \ ( P"b ) ;
% 1 x1=U\ z1 ;
norm( x2-ones ( 3 , 1 ) , i n f ) norm( x1-ones ( 3 , 1 ) , i n f )
% 0 % 0
Uwaga: wyjątkowa sytuacja - brak błędów zaokrągleń. Powtórzyć obliczenia, zmieniając: A2,1 = 1.1, b2 = 6.1.
L.J. Opalski, Wstęp do metod numerycznych - materiały do wykładu wersja: 16.III.2010
Przykład
Rozwiązywanie układów równań liniowych 19
Znalezć transmitancję napięciową K(j�) = U2(j�)/E(j�) liniowej sieci elektrycznej, jak na rysunku, dla wartości
parametrów R1 = 1k&!, R2 = 1M&!, C1 = C2 = 1�F . Wykreślić moduł i fazę transmitancji dla � " [0.1, 105].
R R
1 2
Równania bilansu prądów:
(E - U1)/R1 + (U2 - U1)/R2 = U1 � j�C1
e(t)
u1 u2
C C
(U1 - U2)/R2 = U2 � j�C2
1 2
�ł łł �ł łł �ł łł
1/R1 + 1/R2 + j�C1 -1/R2 �ł �ł U1 �ł �ł E/R1 �ł
�ł
Układ równań liniowych: =
-1/R2 1/R2 + j�C2 U2 0
transmitancja T(j�)
-20
R1=1e3 ; R2=1e6 ; C1=1e-6; C2=C1 ;
-40
AR= [ 1 /R1+1/R2,-1/R2;-1/R2 , 1 / R2 ] ;
-60
AI =[C1, 0 ; 0 , C2 ] ; b = [ 1 /R1 ; 0 ] ;
-80
omega=logspace (-1 ,5);
-100
-120
for n=1: numel ( omega )
-140
A=AR+ j"omega ( n)"AI ;
0 2 4
10 10 10
U=A \ b ;
� [rad/s]
T ( n)=U( 2 ) ;
end
subplot ( 2 , 1 , 1 ) ;
-1
semilogx ( omega,20"log10 ( abs ( T ) ) ) ; axis t i g h t ; grid on ;
ylabel ( |T|[dB] ) ; xlabel ( \omega[rad/s] ) ;
-2
t i t l e ( transmitancjaT(j\omega) ) ;
subplot ( 2 , 1 , 2 ) ; axis t i g h t ; grid on ;
-3
semilogx ( omega , angle ( T ) ) ;
0 2 4
10 10 10
ylabel ( arg(T)[rad] ) ; xlabel ( \omega[rad/s] ) ;
� [rad/s]
L.J. Opalski, Wstęp do metod numerycznych - materiały do wykładu wersja: 16.III.2010
|T| [dB]
arg(T) [rad]
Rozwiązywanie układów równań liniowych 20
4.1.3 Rozkład QR
A = Q � R
R - macierz trójkątna górną
Q - macierz ortogonalna
Ponieważ Q-1 = QT , stąd:
A � x = b
Q � R � x = b
R � x = QT � b
L.J. Opalski, Wstęp do metod numerycznych - materiały do wykładu wersja: 16.III.2010
Rozwiązywanie układów równań liniowych 21
Nakłady obliczeniowe metody QR
rozkład macierzy
2
mnożenie i dzielenie: M H" n3
3
2
dodawanie i odejmowanie: D H" n3
3
pierwiastki kwadratowe: P = n - 1
rozwiązywanie
3
mnożenie i dzielenie: M H" n2
2
3
dodawanie i odejmowanie: D H" n2
2
Rozkład QR w programie MATLAB
1.[Q,R]=qr(A)rozkłada macierzAna iloczynA=QR, przy czymRjest macierzą trójkątną górną, aQ
macierzą unitarną.
2.[Q,R,P]=qr(A)zwraca dodatkowo macierz przestawień kolumnPtaką, żeAP=QR, a elementy wektora
abs(diag(R))są malejące.
L.J. Opalski, Wstęp do metod numerycznych - materiały do wykładu wersja: 16.III.2010
Rozwiązywanie układów równań liniowych 22
Przykład
[Q,R, P]= qr (A)
A=[1 ,1 ,1;1 ,2 ,3;1 , -1 ,1]; % Q =
b = [ 3 ; 1 ; 6 ] ; % -0.3015 -0.2752 -0.9129
[Q,R]= qr (A) % -0.9045 -0.2202 0.3651
% Q = % -0.3015 0.9358 -0.1826
% -0.5774 -0.1543 -0.8018 % R =
% -0.5774 -0.6172 0.5345 % -3.3166 -1.8091 -1.5076
% -0.5774 0.7715 0.2673 % 0 -1.6514 0.4404
% 0 0 -0.7303
% R = % P =
% -1.7321 -1.1547 -2.8868 % 0 0 1
% 0 -2.1602 -1.2344 % 0 1 0
% 0 0 1.0690 % 1 0 0
z=Q "b ; z=Q "b ;
x=R\ z x=P"(R\ z )
% x = % x =
% 1.0000 % 1.0000
% 1.0000 % 1.0000
% 1.0000 % 1.0000
norm( x-ones ( 3 , 1 ) , i n f ) norm( x2-ones ( 3 , 1 ) , i n f )
% 1.7764e-015 % 1.5543e-015
Uwaga: błędy zaokrągleń przy rozwiązywaniu równań dały rozwiązanie o niedokładności 14 EPS.
L.J. Opalski, Wstęp do metod numerycznych - materiały do wykładu wersja: 16.III.2010
Rozwiązywanie układów równań liniowych 23
4.2 Dokładność rozwiązywania układów równań liniowych
4.2.1 Normy wektorów i macierzy
Właściwości normy wektora � w przestrzeni liniowej Rn:
x 0 "x " Rn, x = 0 �! x = 0
ąx = |ą| � x , "ą " R, "x " Rn
x + y x + y , "x, y " Rn
Normy H�ldera (p-normy) wektora z liniowej przestrzeni wektorów N-wymiarowych RN :
1/p
N
x = |xn|p , dla p = 1, 2, . . . , "
p
n=1
N
Ważne p-normy : x = |xn| - norma pierwsza,
1
n=1
N
x = |xn|2 - norma Euklidesa,
2
n=1
x = sup |xn| - norma maksimum (Czebyszewa)
"
n=1,2,...,N
"
Relacje pomiędzy normami: x x x N x N x "
" 2 1 2
L.J. Opalski, Wstęp do metod numerycznych - materiały do wykładu wersja: 16.III.2010
Rozwiązywanie układów równań liniowych 24
Normy macierzy A " RM�N indukowane przez normy wektorów:
Ax p
A = sup = max Ax p
p
x p x =1
x =0
M
A = max |amn|
1
n=1,2,...,N
m=1
Ax 2
A = max = max i(AtA)
2
x =0 x 2 i=1,...,N
gdzie i(B) oznacza wartość własną macierzyB
tzn. liczbę spełniającą równanie Bvi = vi dla pewnego wektora (własnego)vi
N
A " = max |amn|
m=1,2,...,M
n=1
2
Podstawowa relacja pomiędzy normami: A A � A "
1
2
Z definicji normy:
Ax A � x
O normie macierzy i normie wektora, spełniających powyższy warunek, mówimy że są zgodne.
L.J. Opalski, Wstęp do metod numerycznych - materiały do wykładu wersja: 16.III.2010
Rozwiązywanie układów równań liniowych 25
Norma Frobeniusa (Euklidesa) macierzy:
M N
A = |amn|2
E
m=1 n=1
jest zgodna z normą Euklidesa wektorów.
Własności: A A min(M, N) � A 2
2 E
Uwaga:
Realizacja norm macierzy w programieMATLAB:
norma wywołanie sposób realizacji
A1 norm(A,1) max(sum(abs(A)))
A2 norm(A,2) max(svd(A))
AE norm(A, fro ) sqrt(sum(diag(A " A)))
A" norm(A,inf) max(sum(abs(A)))
L.J. Opalski, Wstęp do metod numerycznych - materiały do wykładu wersja: 16.III.2010
Rozwiązywanie układów równań liniowych 26
4.2.2 Niedokładność rozwiązania
" Wpływ zaburzeń wektora b na rozwiązanie x układu równań Ax = b:
A(x + "x) = b + "b �! A � "x = "b �! "x = A-1 � "b
1 1
"x A-1 � "b , oraz b A � x , tzn. A �
x b
Stąd
"x "b
cond(A)
x b
cond(A) = A-1 � A - wskaznik uwarunkowania
" Dla małych zaburzeń, tzn. gdy A-1 � "A < 1, wpływ niedokładności macierzy A można
oszacować następująco:
cond(A) "A
"x
A
x
1 - cond(A) "A
A
Uwaga: wskaznik uwarunkowania liczy w programieMATLABfunkcjacond
L.J. Opalski, Wstęp do metod numerycznych - materiały do wykładu wersja: 16.III.2010
Rozwiązywanie układów równań liniowych 27
Przykład
Dany jest układ równań A � x = b z macierzą Hilberta A, o elementach ai,j = 1/(i + j - 1). Wektor
b dobrano tak, by wektor rozwiązań dokładnych (x) składał się z jedynek.
Propagacja niedokladnosci b Maks. dokl. rozw.
20 5
10 10
||x1-x||
cond(A)
cond(A)*eps
0
T(A) 10
15
10
-5
10
10
10
-10
10
5
10
-15
10
0 -20
10 10
2 4 6 8 10 12 2 4 6 8 10 12
n n
Współczynnik propagacji niedokładności wektora b Niedokładność rozwiązania numerycznego x1, a
szacowano wg: zgrubne oszacowanie - za pomocą: cond(A) " eps.
�
�
x - x b - b
T = /
x b
wykorzystując losowe zaburzenia wektora b (amplituda
�
względna: 0.05%), ozn. b.
L.J. Opalski, Wstęp do metod numerycznych - materiały do wykładu wersja: 16.III.2010
Rozwiązywanie układów równań liniowych 28
4.2.3 Poprawność numeryczna algorytmów skończonych
1. Algorytm Gaussa (oraz metoda LU) z przestawianiem równań i/lub zmiennych są numerycznie
poprawne, a równoważne (błędom zaokrągleń operacji zmiennopozycyjnych) zaburzenie macierzy
układu spełnia warunek:
"A 1
O(n3) � eps
A 1
2. Algorytm rozkładu QR używający przekształceń Hauseholdera jest numerycznie poprawny.
"A E
cn2 � eps, gdzie c H" 4
A E
Dla ułatwienia dowodzenia poprawności algorytmów przydatne jest poniższe
Twierdzenie
Algorytm rozwiązywania układu N równań liniowych jest numerycznie poprawny wtedy i tylko wtedy, gdy
�
istnieje stała K taka, że dla każdego układu równań z tej klasy obliczony wektor reszt: r = fl(Ax - b)
� �
spełnia warunek: r K � A � x � eps , a K = O(N3).
Zadanie (dla zainteresowanych). Pokazać, rozwiązując układy równań z macierzą Hilberta (z
poprzedniego przykładu), że algorytmyMATLABa: LU, QR oraz algorytm obsługujący operator \
zachowują się jak algorytmy poprawne numerycznie.
L.J. Opalski, Wstęp do metod numerycznych - materiały do wykładu wersja: 16.III.2010
Rozwiązywanie układów równań liniowych 29
4.3 Iteracyjne metody rozwiązywania układów równań liniowych
4.3.1 Pochodzenie i własności wielkich zadań
Pochodzenie wielkich układów równań liniowych
" równania dużych sieci (np. elektrycznych)
" dyskretyzacja równań różniczkowych cząstkowych
Cechy wielkich równań liniowych ułatwiające rozwiązywanie:
" małe wypełnienie (procent elementów niezerowych)
" struktura blokowa lub pasmowa macierzy
" dodatnia określoność: "x = 0 : xT Ax > 0

N
" słaba diagonalna dominacja: |amm| |amn|
n=1,n =m
" nieredukowalność macierzy. Macierz kwadratową A nazywamy redukowalną, jeśli istnieje macierz
B C
permutacji P i macierze kwadratowe B, C takie, że P-1AP =
0 C
" znane jest dobre przybliżenie rozwiązania
L.J. Opalski, Wstęp do metod numerycznych - materiały do wykładu wersja: 16.III.2010
Rozwiązywanie układów równań liniowych 30
4.3.2 Wprowadzenie do algorytmów iteracyjnych
Przykład
Należy znalezć rozwiązanie x równania
a � x = 1 , gdzie a jest pewnym parametrem
Równanie to jest równoważne następującemu:
x = 1 + (1 - a) � x
Przyjmijmy jakieś przybliżenie początkowe rozwiązania, np. x(0) = 0;
kolejne przybliżenia 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ązania jest przydatne, jeśli
granica generowanego ciągu jest poszukiwanym rozwiązaniem
w niewielkiej liczbie iteracji można uzyskać dobre przybliżenie rozwiązania.
L.J. Opalski, Wstęp do metod numerycznych - materiały do wykładu wersja: 16.III.2010
Rozwiązywanie układów równań liniowych 31
Zbieżność aperiodyczna Zbieżność oscylacyjna
x(n) dla rownania: 0.5 * x = 1 x(n) dla rownania: 1.5 * x = 1
2 1
0.8
1.5
0.6
1
0 2 4 6 8 10 0 2 4 6 8 10
Numer iteracji (n) Numer iteracji (n)
abs(x(n)-2) abs(x(n)-0.666667)
0 0
10 10
-2 -2
10 10
-4 -4
10 10
0 2 4 6 8 10 0 2 4 6 8 10
Numer iteracji (n) Numer iteracji (n)
Rozbieżność ciągu
x(n) dla rownania: 2.5 * x = 1
50
0
-50
0 2 4 6 8 10
Numer iteracji (n)
abs(x(n)-0.4)
2
10
0
10
-2
10
0 2 4 6 8 10
Numer iteracji (n)
Zachowanie ciągu {x(n)} dla kilku wartości parametru a
L.J. Opalski, Wstęp do metod numerycznych - materiały do wykładu wersja: 16.III.2010
Rozwiązywanie układów równań liniowych 32
4.3.3 Ogólna idea iteracyjnego rozwiązywania układów równań
Zadanie do rozwiązania:
Ax = b
Z macierzy A zostaje wydzielona pewna macierz N
A = N + (A - N)
Nx(i+1) = b - (A - N)x(i)
d(i)
x(i+1) = (I - N-1A) x(i) + N -1b
c
B
Warunki stosowalności:
" układ równań Nx(i+1) = d(i) jest łatwy (tani) do sformułowania i rozwiązania
" ciąg przybliżeń x(i) jest (i to dostatecznie szybko) zbieżny do rozwiązania:
lim x(i) = x" = A-1b
i"
L.J. Opalski, Wstęp do metod numerycznych - materiały do wykładu wersja: 16.III.2010
Rozwiązywanie układów równań liniowych 33
Twierdzenie
Liniowa metoda iteracyjna
x(i+1) = B � x(i) + c
jest zbieżna do rozwiązania wtedy i tylko wtedy, gdy promień spektralny:
�(B) < 1
Ć Ć Ć
oraz zachodzi warunek zgodności: x = Bx + w , gdzie x jest rozwiązaniem dokładnym.
Uwagi:
" �(B) = {max() | " Spect(B)} gdzie widmo Spect macierzy jest zbiorem jej
wartości własnych.
" Szybkość zbieżności jest tym większa im promień spektralny macierzy B jest mniejszy.
L.J. Opalski, Wstęp do metod numerycznych - materiały do wykładu wersja: 16.III.2010
Rozwiązywanie układów równań liniowych 34
4.3.4 Metoda Jacobiego
A = L + U + D
L - macierz dolna trójkątna (z zerową diagonalą)
U - macierz górna trójkątna (z zerową diagonalą)
D - macierz diagonalna (nieosobliwa)
D � x(i+1) = -(L + U) � x(i) + b
d(i)
albo x(i+1) = -D-1(L + U) �x(i) + D-1 � b
c
B
WKiW zbieżności:
" macierz A jest nieredukowalna i diagonalnie słabo dominująca, albo
" macierz A jest silnie diagonalnie dominująca
L.J. Opalski, Wstęp do metod numerycznych - materiały do wykładu wersja: 16.III.2010
Rozwiązywanie układów równań liniowych 35
Przykład
Układ równań do rozwiązania:
�ł łł �ł łł �ł łł
3 -2 0 0 x1 0
�ł śł �ł śł �ł śł
�ł śł �ł śł
-2 5 -2 0 x2 śł �ł 0
�ł śł �ł śł �ł śł
=
�ł śł �ł śł �ł śł
�ł śł �ł śł
0 -2 5 -2 x3 śł �ł 0
�ł �ł �ł �ł �ł �ł
0 0 -2 2 x4 1
x
A=L+D+U b
Równania iteracji metody Jacobiego:
�ł łł �ł łł �ł łł
3 0 0 0 0 -2 0 0 0
�ł śł �ł śł �ł śł
�ł śł �ł śł �ł śł
0 5 0 0 -2 0 -2 0 0
�ł śł �ł śł �ł śł
x(i+1) = - �ł
x(i) +
�ł śł śł �ł śł
�ł śł �ł śł �ł śł
0 0 5 0 0 -2 0 -2 0
�ł �ł �ł �ł �ł �ł
0 0 0 2 0 0 -2 0 1
D L+U b
L.J. Opalski, Wstęp do metod numerycznych - materiały do wykładu wersja: 16.III.2010
Rozwiązywanie układów równań liniowych 36
4.3.5 Metoda Gaussa-Seidla
(L + D) � x(i+1) = -U � x (i) + b
d(i)
albo x(i+1) = -(L + D)(-1) � U �x(i) + (L + D)(-1) � b
B c
WKiW zbieżności:
" macierz A jest dodatnio określona, albo
" macierz A jest silnie diagonalnie dominująca
L.J. Opalski, Wstęp do metod numerycznych - materiały do wykładu wersja: 16.III.2010
Rozwiązywanie układów równań liniowych 37
Przykład
Układ równań do rozwiązania:
�ł łł �ł łł �ł łł
3 -2 0 0 x1 0
�ł śł �ł śł �ł śł
�ł śł �ł śł
-2 5 -2 0 x2 śł �ł 0
�ł śł �ł śł �ł śł
=
�ł śł �ł śł �ł śł
�ł śł �ł śł
0 -2 5 -2 x3 śł �ł 0
�ł �ł �ł �ł �ł �ł
0 0 -2 2 x4 1
A=L+D+U x
b
Równania iteracji metody Gaussa-Seidla:
�ł łł �ł łł �ł łł
3 0 0 0 0 -2 0 0 0
�ł śł �ł śł �ł śł
�ł śł �ł śł �ł śł
-2 5 0 0 0 0 -2 0 0
�ł śł �ł śł �ł śł
x(i+1) = - �ł
x(i) +
�ł śł śł �ł śł
�ł śł �ł śł �ł śł
0 -2 5 0 0 0 0 -2 0
�ł �ł �ł �ł �ł �ł
0 0 -2 2 0 0 0 0 1
L+D U b
L.J. Opalski, Wstęp do metod numerycznych - materiały do wykładu wersja: 16.III.2010
Rozwiązywanie układów równań liniowych 38
4.3.6 Metody relaksacji
1
N(�) = (D + �L)
�
gdzie parametr relaksacji � " (0, 2).
Równania iteracji metod relaksacji:
� - 1
Nx(i+1) = b - (A - N)x(i) = b - (U + � D)x(i)
�
albo x(i+1) = (I - N-1A) x(i) + N -1b
c
B
Metoda kolejnych nadrelaksacji: � > 1
L.J. Opalski, Wstęp do metod numerycznych - materiały do wykładu wersja: 16.III.2010
Rozwiązywanie układów równań liniowych 39
Przykład
Układ równań do rozwiązania:
�ł łł �ł łł �ł łł
3 -2 0 0 x1 0
�ł śł �ł śł �ł śł
�ł śł �ł śł
-2 5 -2 0 x2 śł �ł 0
�ł śł �ł śł �ł śł
=
�ł śł �ł śł �ł śł
�ł śł �ł śł
0 -2 5 -2 x3 śł �ł 0
�ł �ł �ł �ł �ł �ł
0 0 -2 2 x4 1
x
A=L+D+U b
Równania iteracji metod relaksacji:
�ł łł �ł łł �ł łł
3/� 0 0 0 3�-1 -2 0 0 0
�
�ł śł �ł śł �ł śł
�ł śł �ł śł �ł śł
-2 5/� 0 0 0 5�-1 -2 0 0
�ł śł �ł śł �ł śł
�
x(i+1) = - �ł
x(i) +
�ł śł śł �ł śł
�ł śł �ł śł �ł śł
0 -2 5/� 0 0 0 5�-1 -2 0
�ł �ł �ł � �ł �ł �ł
0 0 -2 2/� 0 0 0 2�-1 1
�
1 �-1
b
N=L+ D U+ D
� �
L.J. Opalski, Wstęp do metod numerycznych - materiały do wykładu wersja: 16.III.2010
Rozwiązywanie układów równań liniowych 40
Porównanie szybkości zbieżności algorytmów iteracyjnych
||x(n)-x||
0
10
Jacobi
Gauss-Seidel
SOR, �=1.2
-5
10
-10
10
-15
10
-20
10
0 20 40 60 80 100
Liczba iteracji n
L.J. Opalski, Wstęp do metod numerycznych - materiały do wykładu wersja: 16.III.2010

Wyszukiwarka

Podobne podstrony:
lab 7 zad 0
METROL Zad Domowe2014 15gr29 2 1Mtr
zad domowe redukcja 2
Lab ME II zad rach 12 13
Zdania domowe zad 1 2009 2010
postgreSQL zad lab
W04 zaopatrzenie 2
Załącznik nr 18 zad z pisow wyraz ó i u poziom I
zadanie domowe zestaw
Lab cpp
zad
zestawy domowe ćwiczeń korekcja
lab 2
T2 Skrypt do lab OU Rozdział 6 Wiercenie 3
IE RS lab 9 overview

więcej podobnych podstron