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 nn 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 = = 510-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 zle 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 - zakj
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 wskaznikach 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 wskaznik 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 wskaznikach
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
wskaznika
ę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 lmax 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?
Wyszukiwarka
Podobne podstrony:
Wyklad studport 8Wyklad studport 9Wyklad studport 5Wyklad studport 11Wyklad studport 12Wyklad studport 14Wyklad studport 4Wyklad 12 Podstawowe typy zwiazkow chemicznych blok s i p PCHN SKP studportWyklad 5 Uklad okresowy pierwiastkow studportWyklad 10 Elektrolity, woda, kwasy i zasady PCHN SKP studportWYKŁAD 4 Elektronowa struktura atomu studportSieci komputerowe wyklady dr FurtakWykład 05 Opadanie i fluidyzacjaWYKŁAD 1 Wprowadzenie do biotechnologii farmaceutycznejwięcej podobnych podstron