Oscylacje metody Newtona
Wykres funkcji y = x3 - 3x2 + x + 3
6
5
4
n xn f(xn)
3
0 1.0 2.0
2
1 2.0 1.0
2 1.0 2.0
1
3 2.0 1.0
0
x0
x1
... ... ...
x3
x2
-1
-2
-1 -0.5 0 0.5 1 1.5 2 2.5 3
Rozbieżność metody Newtona
Wykres funkcji y = arc tg (x-2)
1.5
1
n xn f(xn)
0 0.4800 -0.9889
0.5
1 3.7536 1.0525
x2
x0
0
x1
2 -0.5357 -1.1952
x3
3 8.3440 1.4145
-0.5
4 -49.9967 -1.5516
-1
... ... ...
-1.5
-5 0 5 10
Metoda Newtona. Charakterystyka.
" wymaga jednego punktu startowego
" funkcja i jej pierwsza pochodna muszą być ciągłe w przedziale poszukiwań
" warunki zakończenia obliczeń:
" błąd jest dostatecznie mały (długość przedziału zawierającego
pierwiastek jest mniejsza od zadanej dokładności)
" wartość funkcji dostatecznie bliska 0 (jest mniejsza od zadanej
dokładności
" osiągnięto maksymalną zadaną liczbę iteracji
" nie zawsze jest zbieżna
" jest szybsza od metod bisekcji, regula falsi i siecznych (zbieżność kwadratowa)
Metoda Newtona dla układów równań nielinowych. (1)
Analogia dla przypadku
Dwa równania z dwoma niewiadomymi:
pojedynczego równania:
f1(x1, x2) =� 0
równania linearyzujemy, a
f2(x1, x2) =� 0
następnie rozwiązujemy
Zakładamy, że (x1, x2) jest przybliżonym rozwiązaniem układu równań. Obliczamy
poprawki h1i h2 takie, żeby (x1+h1, x2+h2) było lepszym przybliżeniem.
Człony liniowe
ś�f1 ś�f1
0 =� f1(x1 +� h1, x2 +� h2) � f1(x1, x2) +� h1 +� h2
rozwinięcia Taylora
ś�x1 ś�x2
Wszystkie pochodne
ś�f2 ś�f2
0 =� f2(x1 +� h1, x2 +� h2) � f2(x1, x2) +� h1 +� h2
są obliczane w (x1, x2)
ś�x1 ś�x2
ś�f1 ś�f1
�� ł�
ę�ś�x ś�
0 f1(x1, x2) h1 Przypadek jednej zmiennej:
�� ł� �� ł� �� ł�
=� +� ę�ś�f1 ś�x2 ś� ��
ę�0ś� ę�
ś�f2 ś� ę�h2 ś�
ó�
f2(x1, x2)ś� ę� 0 =� f (xn) +� f (xn)(xn+�1 -� xn)
2
�� �� �� �� �� ��
ę� ś�
ś�x2 ��
��ś�x1
Jakobian (oznaczenie J )
Metoda Newtona dla układów równań nielinowych. (2)
Jeżeli macierz J jest nieosobliwa (macierz kwadratowa n��n jest nieosobliwa
jeśli jej rząd jest równy n) to istnieje rozwiązanie:
h1 f1(x1, x2)
�� ł� �� ł�
-�1
=� -�J
ę�h ś� ę�
f2(x1, x2)ś�
2
�� �� �� ��
Zatem
Przypadek jednej zmiennej
( ( (k)
f (xk )
��x1k+�1) ł� ��x1k ) ł� ��h1 ł�
xk +�1 :=� xk -�
:=� +�
ę� ś� ę� ś� ę� ś�
(k+�1) (k ) (k)
ó�
f (xk )
2 2 2
��x �� ��x �� ��h ��
hn
gdzie układ liniowy
(k ( (
��h1 ) ł� ��
f1(x1k ), x2k ) )ł�
J =�
ę� ś� ę� ś�
(k ) ( (
f2(x1k ), x2k ) )��
2
��h �� ��
w praktyce rozwiązujemy metodą eliminacji Gaussa (obliczanie macierzy
odwrotnej jest zbyt kosztowne).
Metoda Newtona dla układów równań nielinowych. (3)
X
f1(x1, x2,K�, xn )
0
�� ł�
�� ł�
ę�
f2(x1, x2,K�, xn )ś� ę�0ś�
ę� ś�
ę� ś�
F(X ) =� 0
=�
ę�M� ś�
ę� ś�
M�
ę� ś�
ę� ś�
fn (x1, x2,K�, xn )�� ��0��
ę� ś�
��
0 =� F(X +� H) � F(X ) +� J(X )H
F
jest odpowiednikiem równania
Wektor poprawek:
ś�f1 ś�f1
-�1 �� ł�
H =� -�J F(X )
ę�ś�x ś�
0 f1(x1, x2) h1
�� ł� �� ł� �� ł�
=� +� ę�ś�f1 ś�x2 ś� ��
W praktyce H wyznaczamy
ę�0ś� ę�
ś�f2 ś� ę�h2 ś�
f2(x1, x2)ś� ę�
2
�� �� �� �� �� ��
metodą eliminacji Gaussa.
ę�
1 ��
��ś�x ś�x2 ś�
Ostatecznie
(k+�1) (k) (k)
X =� X +� H
(k) (k) (k)
gdzie H(k) jest rozwiązaniem układu: J(X )H =� -�F(X )
Metoda Newtona. Układy nielinowe. Przykład. (1)
xy =� z2 +�1 xy -� z2 -�1 =� 0 �� ł�
xy -� z2 -�1
ę� ś�
xyz +� y2 =� x2 +� 2 xyz -� x2 +� y2 -� 2 =� 0
F =�
ę�xyz -� x2 +� y2 -� 2ś�
ę� ś�
x
ex +� z =� ey +� 3 ex -� ey +� z -� 3 =� 0
ę�e -� ey +� z -� 3 ś�
�� ��
y x -� 2z
�� ł�
ę�yz ś�
2
J =� -� 2x xz +� 2x2 xy
�� ł�
x1x2 -� x3 -�1
ę� ś�
ę� ś�
ę� ś�
ex -� ey 1
2 2
�� ��
F =� x2x3 -� x1 +� x2 -� 2ś�
ę�x
1
ę�e -� ex +� x3 -� 3 ś�
x1
x2 x1 -� 2x3 2
�� ł�
�� ��
ę�x
J =� x3 -� 2x1 x1x3 +� 2x2 x1x2 ś�
2
ę� ś�
1 2
ę� ś�
ex -� ex 1
�� ��
function F = uklad3(x)
F = [x(1)*x(2)-x(3)^2-1
x(1)*x(2)*x(3)-x(1)^2+x(2)^2-2
function J=u3jac(x)
exp(x(1))-exp(x(2))+x(3)-3];
J = [ x(2), x(1), -2*x(3)
x(2)*x(3)-2*x(1), x(1)*x(3)+2*x(2), x(1)*x(2)
exp(x(1)), -exp(x(2)), 1 ];
Metoda Newtona. Układy nielinowe. Przykład. (2)
function x = NewtonSys (F, J, x0, tol, kmax)
-�1
H =� -�J F(X )
xold = x0; iter = 1;
while iter <= kmax
(k+�1) (k) (k)
X =� X +� H
h = -feval (J, xold) \ feval (F, xold);
1/ 2
xnew = xold + h';
n
ć� ��
new
��
dif = norm (xnew - xold);
i
��(x -�xiold )2 ��
�� ��
Ł� i=�1 ł�
disp ( [iter xnew dif] );
if dif <= tol
x = xnew;
disp ( 'Metoda Newtona zbiegla sie' )
return
else
x0 = [1 1 1];
xold = xnew;
x = NewtonSys (@uklad3, @u3jac, x0, eps, 10)
end
iter = iter + 1;
end
po 6 iteracjach
disp ( 'Metoda Newtona nie zbiegla sie' )
x = xnew;
x = 1.77767191801074 1.42396059788849 1.23747111773170
Funkcja fzero. (1)
x = fzero (f, x0); -� znajduje zero funkcji f jednej zmiennej w pobliżu punktu
startowego x0
x = fzero (f, x0, opcje); -� znajduje zero funkcji f uwzględniając parametry
podane w strukturze opcje
[x, fval] = fzero (...);
Algorytm  kombinacja
[x, fval, exitflag] = fzero (...);
metod bisekcji, siecznych
[x, fval, exitflag, output] = fzero (...);
oraz odwrotnej interpolacji
kwadratowej
Jeśli x0 jest wektorem dwuelementowym wówczas fzero szuka pierwiastka
w przedziale o końcach określonych przez ten wektor, przy czym funkcja f
musi mieć na końcach przedziału różne znaki.
f jest m-plikiem funkcji albo uchwytem funkcji.
x jest szukanym pierwiastkiem funkcji f.
Funkcja fzero. (2)
Przykład.
-� równanie opisujące głębokość zanurzenia kuli
h3 -� 3h2 +� 4� =� 0
function y = zanurzenie (h);
global ro
global ro
ro = 0.25; % skalar
y = h^3 -� 3*�h^2 + 4*�ro;
x = fzero ('zanurzenie', 0.5);
% wektor
Wykres funkcji y = h3 - 3*h2 + 1
1
x = fzero ('zanurzenie', [0 1]);
0.5
% plik funkcji
uchzan = @zanurzenie;
0
x = fzero (uchzan, 0.5);
% uchwyt funkcji
-0.5
x = fzero (@zanurzenie, 0.5);
-1
0 0.2 0.4 0.6 0.8 1
Funkcja fzero. (3)
[x, fval, exitflag, output] = fzero (f, x0, opcje);
fval  wartość funkcji f w znalezionym punkcie x
exitflag  stan wyjścia z funkcji fzero
1  obliczenia zbieżne do rozwiązania x
 1  obliczenia zostały zakończone przez funkcję wyjścia
 3  wartość funkcji f wynosi NaN albo Inf w przedziale poszukiwań
 4  wartości funkcji f są zespolone w przedziale poszukiwań
 5  obliczenia mogą być zbieżne do punktu osobliwego
output  struktura zawierająca informacje o obliczeniach
opcje  struktura zawierająca parametry dla funkcji fzero
TolX  dokładność zakończenia obliczeń dla zmiennej x
FunValCheck  wyświetla ostrzeżenia gdy funkcja f przyjmuje wartości
zespolone, NaN albo Inf;
'on'  wyświetla
'off' nie wyświetla
Funkcja fzero. (4)
opcje  struktura zawierająca parametry dla funkcji fzero
Display  poziom wyświetlania komunikatów
'off'  nie wyświetla
'iter'  wyświetla dla każdej iteracji
'final'  wyświetla dla końcowej iteracji
'notify'  wyświetla gdy obliczenia nie są zbieżne
Strukturę opcje tworzy się za pomocą funkcji optimset.
Przykład
opcje = optimset ('FunValCheck', 'on', 'TolX', 1e-5);
[x, fval] = fzero (@zanurzenie, 0.5, opcje);
Funkcja fsolve. (1)
x = fsolve (f, x0); -� rozwiązuje układ równań nieliniowych postaci f (x) = 0
gdzie x jest wektorem, a f  funkcją wektorową. Punkt x0
jest punktem startowym.
x = fsolve (f, x0, opcje); -� uwzględnia parametry podane w strukturze opcje
[x, fval] = fsolve (...);
Algorytm  korzysta z
[x, fval, exitflag] = fsolve (...);
metod Gaussa-Newtona,
[x, fval, exitflag, output] = fsolve (...);
Levenberga-Marquardta
oraz metody large-scale,
[x, fval, exitflag, output, jakobian] = fsolve (...);
opartych na nieliniowej
metodzie najmniejszych
f  m-plik funkcji albo uchwyt funkcji
kwadratów
x  szukany pierwiastek funkcji f
fval  wartość funkcji f w znalezionym punkcie x
exitflag  stan wyjścia z funkcji fzero
output  struktura zawierająca informacje o obliczeniach
jakobian  macierz Jakobiego w znalezionym punkcie x
Funkcja fsolve. Przykład 1
Rozwiązać układ równań
1
��
2x1 -� x2 -� e-�x =� 0
1
��
��
2x1 -� x2 =� e-�x
2
��-� x1 +� 2x2 -� e-�x =� 0
��
2
��-� x1 +� 2x2 =� e-�x
Punkt startowy
function F = nieliniowy (x);
x0 = [-5; -5];
F = [ 2*x(1) - x(2) - exp(-x(1));
[x, fval] = fsolve (@nieliniowy, x0);
-x(1) + 2*x(2) - exp(-x(2)) ];
x =
0.5671
Jacobian liczony numerycznie
0.5671
fval =
1.0e-006 *
-0.4059
-0.4059
Funkcja fsolve. Przykład 2 (a)
xy =� z2 +�1 xy -� z2 -�1 =� 0 �� ł�
xy -� z2 -�1
ę� ś�
xyz +� y2 =� x2 +� 2 xyz -� x2 +� y2 -� 2 =� 0
F =�
ę�xyz -� x2 +� y2 -� 2ś�
ę� ś�
x
ex +� z =� ey +� 3 ex -� ey +� z -� 3 =� 0
ę�e -� ey +� z -� 3 ś�
�� ��
x2 x1 -� 2x3
�� ł�
ę�x
J =� x3 -� 2x1 x1x3 +� 2x2 x1x2 ś�
2
ę� ś�
1 2
ę� ś�
ex -� ex 1
function [F, J] = uklad3jac(x)
�� ��
F = [x(1)*x(2)-x(3)^2-1
x(1)*x(2)*x(3)-x(1)^2+x(2)^2-2
Jacobian liczony
exp(x(1))-exp(x(2))+x(3)-3];
analitycznie
J = [ x(2), x(1), -2*x(3)
x(2)*x(3)-2*x(1), x(1)*x(3)+2*x(2), x(1)*x(2)
exp(x(1)), -exp(x(2)), 1 ];
Funkcja fsolve. Przykład 2 (b)
opcje = optimset ('Jacobian', 'on', 'Display', 'iter', 'NonlEqnAlgorithm', 'gn');
[x, Fval]=fsolve (@uklad3, [1; 1; 1], opcje)
Iteration Func-count Residual Step-size derivative
Algorytm
0 1 6
Gaussa-Newtona
1 5 0.116285 0.592 -2.81
2 10 3.49932e-006 0.97 -2e-005
3 14 3.13957e-014 1 -8.14e-015
4 15 5.85454e-023 1 -6.28e-014
Optimization terminated: directional derivative along
search direction less than TolFun and infinity-norm of
gradient less than 10*(TolFun+TolX).
x = 1.77767187811036
Otrzymane za pomocą
1.42396056622665
skryptu NewtonSys:
1.23747112327120
Fval = 1.0e-011 * x = 1.77767191801074
0.54756199574513 1.42396059788849
0.32058800059076 1.23747111773170
0.42761350016463
Równowaga chemiczna. Przykład 3 (a)
Znajdz
równowagowe stopnie przereagowania dla układu następujących reakcji:
CC
Początkowe stężenia:
2A +� B �� C K1 =� =� 5��10-�4
2
CACB
CA0 = 40 kmol/m3
CC
A +� D �� C K2 =� =� 4 ��10-�2
CB0 = 15 kmol/m3
CACD
CC0 = 0 kmol/m3
Oznaczenia:
CD0 = 10 kmol/m3
x1  stopień przereagowania składnika B
x2  stopień przereagowania składnika D
Układ równań:
Stężenia składników zależą od stopni konwersji:
CC
CA =� CA0 -� 2x1CB0 -� x2CD0
f1(x1, x2) =� -� K1
2
CACB
CB =� (1-� x1)CB0
CC
CC =� CC0 +� x1CB0 +� x2CD0 f2(x1, x2) =� -� K2
CACD
CD =� (1-� x2)CD0
Równowaga chemiczna. Przykład 3 (b)
function f = reakcja (x, c0, K)
cA = c0(1) - 2*x(1)*c0(2) - x(2)*c0(4);
cB = (1 - x(1))*c0(2);
cC = c0(3) + x(1)*c0(2) + x(2)*c0(4);
Wyniki:
cD = (1 - x(2))*c0(4);
x =
f(1) = cC/cA^2/cB - K(1);
0.1203 0.4787
f(2) = cC/cA/cD - K(2);
f = f ';
fval =
1.0e-006 *
-0.0001
k = [5e-4 4e-2];
0.7938
c0 = [40 15 0 10];
x0 = [0.1 0.9];
[x, fval] = fsolve (@ (x) reakcja (x, c0, k), x0)
Układy równań liniowych. Wstęp.
a11x1 +� a12x2 +�K�+� a1nxn =� b1
a21x1 +� a22x2 +�K�+� a2nxn =� b2
M� M�
an1x1 +� an2x2 +�K�+� annxn =� bn
a11 a12 K� a1n x1 b1
�� ł� �� ł� �� ł�
ę�a
a22 K� a2n ś� ę�x2 ś� ę�b2 ś�
21
ę� ś� ę� ś� ę� ś�
=�
ę� ś� ę� ś� ę� ś�
M� M� M� M� M�
ę� ś� ę� ś� ę� ś�
an2 K� ann�� ��xn �� ��bn ��
n1
��a
b
x
A
Ax =� b
System mieszalników (1)
Pierwszy mieszalnik:
F01c01 +� F31c3 =� F12c1 +� F15c1
F55 = 2
Podstawiając:
c5
F15 = 3
6c1 -� c3 =� 50
F54 = 2
F25 = 1
F12 = 3
F01 = 5 F24 = 1 F44 = 11
c1 c2 c4
c01 = 10
F23 = 1
F31 = 1 F34 = 8
c3
F03 = 8
c  stężenie w mg/m3
c03 = 20
F  natężenie przepływu w m3/min
System mieszalników (2)
11.51
6c1 -� c3 =� 50
�� ł�
ę�11.51ś�
-� 3c1 +� 3c2 =� 0
ę� ś�
ę� ś�
c =� 19.06
-� c2 +� 9c3 =�160
ę� ś�
-� c2 -� 8c3 +�11c4 -� 2c5 =� 0
ę�17.00ś�
ę�11.51ś�
-� 3c1 -� c2 +� 4c5 =� 0
�� ��
Układy dwóch równań. Własności rozwiązań.
ć� ��
a11 �� b1
��
x2 =� -��� +�
a12 ��x1 a12
a11x1 +� a12x2 =� b1
Ł� ł�
a21x1 +� a22x2 =� b2
ć� ��
a21 �� b2
��
x2 =� -��� +�
a22 ��x1 a22
Ł� ł�
Można przedstawić jako proste we
współrzędnych kartezjańskich
a11x1 +� a12x2 =� b1
a21x1 +� a22x2 =� b2
Przykład. Układ nieosobliwy.
3x1 +� 2x2 =�18 9
8
-� x1 +� 2x2 =� 2
7
3x1 +� 2x2 =�18
6
5
4
3
-� x1 +� 2x2 =� 2
2
1
0
0 1 2 3 4 5 6
x1
2
x
Przykład. Układ osobliwy. Brak rozwiązania
1
-� x1 +� x2 =�1
2 2.5
1 1
-� x1 +� x2 =�
2
2 2
1.5
1
-� x1 +� x2 =�1
2
1
0.5
1 1
-� x1 +� x2 =�
0
2 2
-0.5
-1
0 1 2 3 4 5 6
x1
2
x
Przykład. Układ osobliwy. Nieskończenie wiele rozwiązań
1
-� x1 +� x2 =�1
2.5
2
-� x1 +� 2x2 =� 2
2
1.5
1
-� x1 +� x2 =�1
1 2
0.5
-� x1 +� 2x2 =� 2
0
-0.5
-1
0 1 2 3 4 5 6
x1
2
x
Przykład. Układ z
le uwarunkowany
1
-� x1 +� x2 =�1
26.5
2
-� 0.46x1 +� x2 =�1.1
26
-� 0.46x1 +� x2 =�1.1
25.5
25
24.5
1
-� x1 +� x2 =�1
2
24
50 51 52 53 54 55
x1
2
x
Układy łatwe do rozwiązania (1)
a11 0 K� 0 x1 b1 x1 b1 / a11
�� ł� �� ł� �� ł� �� ł� �� ł�
ę� ś� ę�x ś� ę�b ś� ę�x ś� ę�b / a22ś�
0 a22 K� 0
2 2 2 2
ę� ś� ę� ś� ę� ś� ę� ś� ę� ś�
=� =�
ę� ś� ę� ś� ę� ś� ę� ś� ę� ś�
M� M� M� M� M� M� M�
ę� ś� ę� ś� ę� ś� ę� ś� ę� ś�
0 0 K� ann�� ��xn �� ��bn ��
�� n n
��x �� ��b / ann��
Jeżeli dla pewnego i jest aii = 0 oraz bi = 0 to xi może być dowolną liczbą.
Jeżeli natomiast aii = 0 oraz bi ą� 0 to układ nie ma rozwiązania.
Układy łatwe do rozwiązania (2)
a11x1 =� b1
a11 0 K� 0 x1 b1
�� ł� �� ł� �� ł�
ę�a ś� ę�x ś� ę�b ś�
a21x1 +� a22x2 =� b2
a22 K� 0
21 2 2
ę� ś� ę� ś� ę� ś�
=�
ę� ś� ę� ś� ę� ś� M�
M� M� M� M� M�
ę� ś� ę� ś� ę� ś�
an1x1 +� an2x2 +�K�+� annxn =� bn
an2 K� ann�� ��xn �� ��bn ��
n1
��a
Algorytm podstawienia w przód:
x1 =� b1 / a11
input n, (aij), (bi)
x2 =� (b2 -� a21x1) / a22
for i = 1 to n do
s Ź� bi
M�
for j = 1 to i-1 do
xn =� (bn -� an1x1 -� an2x2-�K�-� an(n-�1)xn-�1) / ann
s Ź� s -� aij xj
end do
xi Ź� s / aii
i-�1
end do
xi Ź� (bi -�
output (xi)
ij
��a xj ) / aii
j=�1
Układy łatwe do rozwiązania (3)
a11 a12 K� a1n x1 b1
�� ł� �� ł� �� ł�
a11x1 +� a12x2 +�K�+� a1nxn =� b1
ę�
0 a22 K� a2n ś� ę�x2 ś� ę�b2 ś�
a22x2 +�K�+� a2nxn =� b2
ę� ś� ę� ś� ę� ś�
=�
ę� ś� ę� ś� ę� ś�
M� M� M� M� M�
M� M�
ę� ś� ę� ś� ę� ś�
0 0 K� ann�� ��xn �� ��bn ��
annxn =� bn
��
xn =� bn / ann
Algorytm podstawienia wstecz:
xn-�1 =� (bn-�1 -� a(n-�1)nxn ) / a(n-�1)(n-�1)
input n, (aij), (bi)
for i = n to 1 step -1 do
M�
s Ź� bi
x1 =� (b1 -� a1nxn -� a1(n-�1)xn-�1-�K�-� a12x2) / a11
for j = i+1 to n do
s Ź� s -� aij xj
end do
n
xi Ź� s / aii
xi Ź� (bi -�
end do
ij
��a xj ) / aii
j=�i+�1
output (xi)
Układy łatwe do rozwiązania (4)
Załóżmy, że znamy permutację
1 a11 a12 0 x1 b1
�� ł� �� ł� �� ł�
(p1, p2, ..., pn) układu (1, 2, ..., n).
ę�a ę�b ś�
2 a22 a23ś� ę�x2 ś� =�
21 2
ę� ś� ę� ś� ę� ś�
Zmodyfikowany algorytm
ę� ś� ę� ś� ę� ś�
3 0 0
31 3 3 podstawienia w przód:
��a �� ��x �� ��b ��
Przestawiamy input n, (aij), (bi), (pi)
wiersze for i = 1 to n do
s Ź� bpi
3 a31 0 0 x1 b3
�� ł� �� ł� �� ł�
for j = 1 to i-�1 do
ę�a ś� ę�x ś� ę�b ś�
s Ź� s -� api j xj
1 a12 0 =�
11 2 1
ę� ś� ę� ś� ę� ś�
end do
ę� ś� ę� ś� ę� ś�
2 a22 a23�� ��x3 �� ��b2 ��
21 xi Ź� s / api i
��a
end do
output (xi)
i-�1
xi Ź� (bp -�
pi j i
��a xj ) / ap
i i
Podobnie można zmodyfikować
j=�1
algorytm podstawienia wstecz.
Rozkłady LU
L  macierz
Ax = b
trójkątna dolna
Jeżeli macierz A można przedstawić w formie A = LU
U  macierz
LUx = b
trójkątna górna
z Ź� wektor
wtedy rozwiązanie dzieli się
na dwa łatwe etapy:
Lz = b oraz Ux = z
Jak to zrobić ?
Eliminacja Gaussa (1)
4 12 8 4
�� ł�
ę�1 7 18 9 ś�
ę� ś�
A =�
ę� ś�
2 9 20 20
ę� ś�
��3 11 15 14��
Wartości początkowe:
4 12 8 4
1 0 0 0 �� ł�
�� ł�
ę�1 7 18 9 ś�
ę�0 1 0 0ś�
ę� ś�
ę� ś�
U =�
L =�
ę� ś�
ę� ś�
2 9 20 20
0 0 1 0
Odejmujemy 1-szy wiersz
ę� ś�
ę� ś�
pomnożony przez:
��3 11 15 14��
��0 0 0 1��
" a21/a11 od 2-go wiersza
Krok 1:
" a31/a11 od 3-go wiersza
4 12 8 4
1 0 0 0 �� ł�
�� ł�
ę�0 4 16 8 ś�
ę�1/ 4 1 0 0ś�
" a41/a11 od 4-go wiersza
ę� ś�
ę� ś�
U =�
L =�
ę� ś�
ę� ś�
0 3 16 18
1/ 2 0 1 0
ę� ś�
ę� ś�
��0 2 9 11��
��3/ 4 0 0 1��
Eliminacja Gaussa (2)
1 0 0 0 4 12 8 4
�� ł� �� ł�
ę�1/ 4 1 0 0ś� ę�0 4 16 8 ś�
ę� ś� ę� ś�
L =� U =�
ę� ś� ę� ś�
1/ 2 0 1 0 0 3 16 18
ę� ś� ę� ś�
��3/ 4 0 0 1�� ��0 2 9 11��
Odejmujemy 2-gi wiersz
pomnożony przez:
Krok 2:
" a32/a22 od 3-go wiersza
4 12 8 4
1 0 0 0 �� ł�
�� ł�
ę�0 4 16 8 ś�
ę�1/ 4 1 0 0ś�
" a42/a22 od 4-go wiersza
ę� ś�
ę� ś�
U =�
L =�
ę� ś�
ę� ś�
0 0 4 12
1/ 2 3/ 4 1 0
ę� ś�
ę� ś�
��0 0 1 7 ��
��3/ 4 1/ 2 0 1��
Odejmujemy 3-ci wiersz
Krok 3:
pomnożony przez:
4 12 8 4
1 0 0 0 �� ł�
�� ł�
" a43/a33 od 4-go wiersza
ę�0 4 16 8 ś�
ę�1/ 4 1 0 0ś�
ę� ś�
ę� ś�
U =�
L =�
ę� ś�
ę� ś�
0 0 4 12
1/ 2 3/ 4 1 0
ę� ś�
ę� ś�
��0 0 0 4 ��
��3/ 4 1/ 2 1/ 4 1��
Eliminacja Gaussa (3)
4 12 8 4
4 12 8 4 �� ł�
�� ł�
1 0 0 0
�� ł�
ę�1 7 18 9 ś�
ę�0 4 16 8 ś�
ę�1/ 4 1 0 0ś�
ę� ś�
ę� ś�
ę� ś�
A =�
U =�
L =�
ę� ś�
ę� ś�
2 9 20 20
0 0 4 12
ę� ś�
1/ 2 3/ 4 1 0
ę� ś�
ę� ś�
ę� ś�
��3 11 15 14��
��0 0 0 4 ��
��3/ 4 1/ 2 1/ 4 1��
input n, (aij)
Twierdzenie:
for k = 1 to n-�1 do
Jeżeli wszystkie elementy główne akk for i = k +1 to n do
obliczane w opisany wyżej sposób są z Ź� aik / akk
różne od 0, to A = LU. aik Ź� 0
for j = k +1 to n do
aij Ź� aij -� z��akj
end do
end do
output (ai j)
Eliminacja Gaussa. Trudności. (1)
Przykład 1. Przykład 2.
��1 e� -�1ł� x1 ��e� -�1ł�
�� ł�
0 1 x1 1
�� ł� �� ł� �� ł�
=�
ę� ś� ę� ś�
ę�x ś�
=�
ę�1 1ś� ę�x ś� ę�2ś�
2
2
�� ��
��1 1 �� �� ��
2
�� �� �� �� �� ��
-�1
e� jest małą liczbą, więc e� jest
Algorytm Gaussa zawodzi
dużą liczbą. Z algorytmu Gaussa
gdyż a11 = 0.
otrzymujemy
Rozwiązanie istnieje -�1
��1 e� -�1 ł� x1 �� ł�
�� ł�
e�
=�
x2 = 1 oraz x1 = 1.
ę� ś� ę� ś�
ę�x ś�
-�1 -�1
2
�� ��
��0 1-� e� �� ��2 -� e� ��
Rozwiązanie
-�1
2 -� e�
-�1
Dokładne rozwiązanie x2 =� x1 =� (1-� x2)e�
-�1
1-� e�
1 1-� 2e�
W arytmetyce komputerowej
x1 =� �1 x2 =� �1
-�1
1-� e� 1-� e� licznik i mianownik są równe -�e�
czyli x2 = 1 oraz x1 = 0.
Eliminacja Gaussa. Trudności. (2)
Przykład 3.
Przestawienie wierszy
e� 1 x1 1
�� ł� �� ł� �� ł�
1 1 x1 2
�� ł� �� ł� �� ł�
=�
ę�1 1ś� ę�x ś� ę�2ś�
=�
ę�e� 1ś� ę�x ś� ę�1ś�
2
�� �� �� �� �� ��
�� �� �� 2 �� �� ��
Stosując algorytm Gaussa
e� jest małą liczbą różną od 0.
dostajemy
Algorytm Gaussa daje
1 1 x1 2
�� ł� �� ł� �� ł�
e� 1 x1 1 =�
�� ł� �� ł� �� ł�
ę�0 1-�e� ś� ę�x ś� ę�1-� 2e� ś�
=�
ę�0 1-� e� ś� ę�x ś� ę�2 -� e� ś�
-�1 -�1
�� �� �� 2 �� �� ��
�� �� �� 2 �� �� ��
Rozwiązania komputerowe
Rozwiązanie
x2 = 1 oraz x1 = 2  x2 = 1.
-�1
2 -� e�
-�1
x2 =� x1 =� (1-� x2)e�
-�1
1-� e�
Dokładne rozwiązanie
W arytmetyce komputerowej
-�1
licznik i mianownik są równe -�e�
1 1-� 2e�
x1 =� �1 x2 =� �1
czyli x2 = 1 oraz x1 = 0.
1-� e� 1-� e�
Algorytm z przestawianiem wierszy
Algorytm musi uwzględniać przestawianie wierszy, jeżeli wymagają tego
okoliczności. Przestawianie wierszy w pamięci komputera wydłuża czas
obliczeń. Aby tego uniknąć inaczej je wybieramy. Zamiast w porządku
naturalnym 1, 2, ..., n-�1, bierzemy kolejno wiersze o wskaz
nikach p1, p2, ..., pn-1.
input n, (ai j), (pi)
for k = 1 to n-�1 do
Algorytm różni się od
for i = k +1 to n do
poprzedniego tylko tym, że
z Ź� api k / api k
pierwszy wskaz
nik tablicy a
api k Ź� 0
zależy od wyboru permutacji.
for j = k +1 to n do
api j Ź� api j -� z apk j
Jak wybrać permutację?
end do
end do
output (ai j)
Skalowany wybór wierszy głównych (1)
k = 1
a11 / s1 =� 2/ 6 2 3 -� 6 6 1
�� ł� �� ł� �� ł�
gdzie:
ę�1 -� 6 8 ś� ę�8ś� ę�2ś�
A =� s =� p =�
a21 / s2 =�1/8
ę� ś� ę� ś� ę� ś�
si =� maxaij
ę� ę� ę� 1Ł� jŁ�n
��3 -� 2 1 ś� ��3ś� ��3ś�
�� �� ��
a31 / s3 =� 3/ 3
Wyznaczenie
największy
3 -� 2 1 3 3
�� ł� �� ł� �� ł�
skali dla
ę�1 -� 6 8 ś� ę�8ś� ę�2ś�
A =� s =� p =�
każdego
ę� ś� ę� ś� ę� ś�
ę� ś� ę� ę� wiersza
��2 3 -� 6�� ��6ś� ��1ś�
�� ��
Odejmujemy 1-szy wiersz
��3 -� 2 1 ł�
pomnożony przez:
ę�0 -� 16 23 ś�
" a21/a11 od 2-go wiersza
3 3
ę� ś�
ę�0 13 -� 20ś�
" a31/a11 od 3-go wiersza
3 3
�� ��
Skalowany wybór wierszy głównych (2)
k = 2
3 3
�� ł� �� ł�
��3 -� 2 1 ł�
ę�8ś� ę�2ś�
ę�0 -� 16 23 ś�
16 2
s =� p =�
a22 / s2 =� /8 =�
3 3
3 3
ę� ś� ę� ś� ę� ś�
13 13
ę�0 13 -� 20ś�
ę� ę�
a32 / s3 =� / 6 =�
�� ��
3 18 ��6ś� ��1ś�
3 3
�� ��
3
�� ł�
��3 -� 2 1 ł�
największy
ę�1ś�
ę�0 36 -� 20ś�
p =�
3 3
ę� ś� ę� ś�
ę�0 -� 16 23 ś�
ę�
��2ś�
��
3 3
�� ��
permutacja
��3 -� 2 1 ł�
ę�0 36 -� 20ś�
U =�
Odejmujemy 2-gi wiersz
3 3
ę� ś�
7 pomnożony przez:
ę�0 0 -� 13ś�
�� ��
" a32/a22 od 3-go wiersza
Skalowany wybór wierszy głównych (3)
Obliczenia dają nam także macierz permutacji:
Macierz P powstaje z
3 0 0 1
�� ł� �� ł�
macierzy jednostkowej I
ę�1 0 0ś�
ę�1ś�
P =�
p =�
ę� ś�
ę� ś�
przez przestawienie wierszy
ę� ś�
ę�
��0 1 0��
��2ś�
��
zgodnie z tablicą p.
oraz rozkład LU macierzy PA:
0 0 1 2 3 -� 6 3 -� 2 1 ��1 0 0ł� ��3 -� 2 1 ł�
�� ł� �� ł� �� ł�
ę�2 1 0ś� ę�0 13 -� 20ś�
ę�1 0 0ś� ę�1 -� 6 8 ś� ę�2 3 -� 6ś�
PA =� =� =�
3 3 3
ę� ś� ę� ś�
ę� ś� ę� ś� ę� ś�
16 7
ś�
ę� ś�
1ś� ę�0 0 -�
��0 1 0�� ę� -� 2 1 ś� ę� -� 6 8 ś� ę�1 -�
��3 �� ��1 ��
3
�� 13 �� �� 13 ��
PA
P A L U
Zamiast rozkładu LU macierzy A otrzymaliśmy rozkład LU macierzy PA.
Skalowany wybór wierszy głównych (4)
W rzeczywistości nie przestawiamy wierszy w pamięci komputera, ale zamiast
w porządku naturalnym 1, 2, ..., n-�1, bierzemy kolejno wiersze o wskaz
nikach
p1, p2, ..., pn-1.
2
13 20
��0 -� ł� 3
�� ł�
3 3
2 3 -� 6
�� ł�
1
�� ł�
ę�0 0 -� 7 ś�
ę�1ś�
ę�1 -� 6 8 ś�
ę�2ś�
A =�
p =�
A =�
p =�
13
ę� ś�
ę� ś� ę� ś�
ę� ś�
ę�
ę�3 2 ś�
ę�
1
ę�
��3 -� 2 1 ś�
��
��3ś�
��
�� ��
��2ś�
��
3
1
13 20
��0 -� ł�
3
�� ł�
3 3
ę�0 -� 16 23 ś�
ę�2ś�
A =�
p =�
3 3
ę� ś�
ę� ś�
ę�3 2 ś�
1
ę�
�� ��
��1ś�
��
Skalowany wybór wierszy głównych (5)
Dodatkowo zapamiętujemy mnożniki na miejscu zerowanych elementów
macierzy A. Zmieniona macierz zawiera wszystkie wielkości potrzebne do
odtworzenia rozkładu LU. Niszczymy jednak pierwotną macierz.
2 13 20
��[� ]� -� ł� 3
�� ł�
3 3 3
2 3 -� 6
�� ł�
1
�� ł�
ę� ś�
ę�1ś�
1 16 7
ę�1 -� 6 8 ś�
ę�2ś�
A =� [� ]� [�-� ]� -�
p =�
A =�
p =�
3 13 13
ę� ś�
ę� ś� ę� ś�
ę� ś�
ę� ę� ś�
ę�
3 2 1
ę�
��3 -� 2 1 ś�
��
��3ś�
��
�� ��
��2ś�
��
2 13 20
��[� ]� -� ł�
3
�� ł�
3 3 3
ę� ś�
1 16 23 ę�2ś�
A =� [� ]� -�
p =�
3 3 3
ę� ś�
ę� ś�
ę� ś�
3 2 1
ę�
�� ��
��1ś�
��
Skalowany wybór wierszy głównych (6)
Algorytm rozwiązywania układu równań Ax = b metodą eliminacji Gaussa
ze skalowanym wyborem wierszy głównych. Etapy:
1) Znalezienie permutacji (p1, p2, ..., pn) zbioru (1, 2, ..., n) oraz rozkładu
LU macierzy PA.
2) Rozwiązanie dwóch układów równań Lz = Pb (dostajemy wektor z) oraz
Ux = z (otrzymujemy wektor x).
Ax = b
PAx = Pb, PA = LU
LUx = Pb, Ux = z (oznaczamy)
Lz = Pb, (algorytm podstawienia w przód)
Ux = z, (algorytm podstawienie wstecz)
Algorytm rozkładu LU
input n, (ai j)
for i = 1 to n do
pi Ź� i (początkowa naturalna kolejność)
si Ź� maxaij
(wyznaczenie skali dla wierszy)
1Ł� jŁ�n
end do
for k = 1 to n-�1 do
wybór takiego j ł� k, że |apj k| / spj ł� |api k| / spi dla i = k, k+1, ..., n
pk Ź� pj
for i = k +1 to n do
z Ź� api k / apk k (wyznaczamy mnożnik z dla wiersza pi)
api k Ź� z (zapamiętujemy mnożnik z)
for j = k +1 to n do
api j Ź� api j -� z apk j (odejmujemy wiersz pk od pi)
end do
end do
end do
output (ai j), (pi)
Algorytm rozwiązania układu Ax = b
po znalezieniu rozkładu PA = LU
input n, (aij), (pi), (bi)
for k = 1 to n-�1 do
for i = k+1 to n do
bpi Ź� bpi -� api k bpk (Lz = Pb, podstawienie w przód)
end do
end do
for i = n to 1 step -�1 do
xi Ź� bpi (Ux = z, podstawienie wstecz)
n
for j = i+1 to n do
xi Ź� (bp -�
xi Ź� xi -� api j xj
pi j i
��a xj ) / ap
i i
j=�i+�1
end do
end do
xi Ź� xi / api i
end do
output (xi)
Koszt obliczeń (1)
dodawanie i odejmowanie  trwają krócej niż mnożenie i dzielenie
mnożenie i dzielenie  zajmują podobną ilość czasu
działania długie: para mnożenie  dodawanie
Koszt obliczeń jest mierzony ilością długich działań.
Pierwszy wiersz:
Rozkład macierzy n-tego stopnia:
obliczenie mnożnika
Krok 1:
-� 1 działanie
określenie a11 a12 K� a1n
�� ł�
a11 / s1
odjęcie wielokrotności
wskaz
nika
ę�a
wiersza 1-go od i-tego
a22 K� a2n ś�
a21 / s2
21
p1 wiersza
ę� ś�
-� (n -� 1) działań
głównego: ę� ś�
M� M� M�
M�
łącznie n operacji
ę� ś�
n działań
an1 / sn
an2 K� ann��
n1
��a
Dla (n -� 1) wierszy:
n(n -� 1) działań
Razem dla 1-go kroku: n + n(n -� 1) = n2 działań
Koszt obliczeń (2)
Dla 1-go kroku: n2 działań
Dla 2-go kroku: (n -� 1) 2 działań
...
Dla n-�1-go kroku: 22 działań
Dla całego rozkładu:
1 1 1 1 1
n2 +� (n -�1)2 +�K�+� 32 +� 22 =� n3 +� n2 +� n � n3 +� n2
3 2 6 3 2
1 1
Rozwiązanie układu Lz = Pb: n2 -� n
2 2
1 1
n2 +� n
Rozwiązanie układu Ux = z:
2 2
Razem:
Obliczenie A-�1 jest
Rozwiązanie układu Ax = b
kosztowne. Wymaga
metodą eliminacji Gaussa ze
1 3
n3 +� n2
4
skalowanym wyborem
3 2
Oć� n3 �� działań
�� ��
3
elementów głównych:
Ł� ł�
Normy wektorów
Każdemu wektorowi x można przyporządkować liczbę ||x|| nazywana normą
wektora x. Liczby te muszą spełniać następujące aksjomaty:
|| x || ł� 0 dla każdego x; || x || = 0 wtedy i tylko wtedy gdy x = 0
|| l� x || = |l�| || x || dla wszystkich x oraz wszystkich l� �� ��
|| x + y || Ł� || x || + || y || dla wszystkich x i y (nierówność trójkąta)
Suma wartości bezwzględnych
Euklidesowa: Maksimum:
współrzędnych:
1/ 2
n
n
ć� ��
2
��
x =�
x =� xi
x =� max xi
i
��x ��
��
2
�� �� 1
Ą�
1Ł�iŁ�n
i=�1
Ł� ł�
i=�1
|| x ||2 Ł� 1 || x ||Ą� Ł� 1 || x ||1 Ł� 1
||
x
-�
y
||
-�
odległość
Normy macierzy
Normy macierzy definiuje się tak samo jak dla wektorów tzn. muszą one
spełniać te same aksjomaty.
Celowy jest taki wybór normy macierzowej || A ||, aby dla danej normy
wektorowej || x || spełniona była nierówność
|| Ax || Ł� || A || || x ||
Takie normy nazywamy zgodnymi.
1/ 2
n
gdzie l�max jest największą
ć� ��
2
��
x =�
A =� max
wartością własną macierzy
i
��x ��
2
�� ��
2
Ł� i=�1 ł�
ATA
(norma spektralna)
n
x =� max xi
(norma wierszowa)
A =� max aij
Ą�
��
Ą�
1Ł�iŁ�n 1Ł�iŁ�n
j=�1
n
n
(norma kolumnowa)
A =� max aij
x =� xi
��
�� 1
1
1Ł� jŁ�n
i=�1
i=�1
Do czego przydają się normy?