Lucrarea de laborator nr. 5 · 2019. 10. 6. · Mădălina Roxana Buneci Metode Numerice...
Transcript of Lucrarea de laborator nr. 5 · 2019. 10. 6. · Mădălina Roxana Buneci Metode Numerice...
-
Metode Numerice - Lucrarea de laborator 5
1
Lucrarea de laborator nr. 5
I. Scopul lucrării
Rezolvarea ecuațiilor/sistemelor neliniare
II. Conținutul lucrării
1. Familia de comenzi solve din MAPLE
2. Metoda bisecției
3. Metoda coardei
4. Metoda tangentei
5. Metoda punctului fix
6. Metoda Newton – cazul m-dimensional
III. Prezentarea lucrării
III.1. Familia de comenzi solve din MAPLE
Spre deosebire de atribuiri, ecuațiile sunt expresii matematice simple care
stabilesc relații între anumite variabile și/sau valori (fără a asocia vreo valoare explicită
pentru variabilele conținute). Operatorul folosit este = (reamintim că în cazul atribuirii
se utilizează : =). O familie de comenzi care folosesc ecuațiile ca parametrii este
familia de comenzi solve. Forma generală a comenzii solve este:
>solve(ecuatie, necunoscuta);
sau
>solve(ecuatii, neconoscute);
MAPLE dispune de comenzi specializate pentru rezolvarea diverselor tipuri de ecuații
(rezolvarea ecuațiilor în diverse mulțimi):
fsolve: rezolvă ecuații aplicând aritmetica virgulei mobile
isolve: rezolvă ecuații în mulțimea numerelor întregi
msolve: rezolvă ecuații modulo m
-
Mădălina Roxana Buneci Metode Numerice –Laborator
2
rsolve: rezolvă ecuații date prin relații de recurență
dsolve: rezolvă ecuații diferențiale ordinare
pdsolve: rezolvă ecuații cu derivate parțiale
linsolve: rezolvă sisteme de ecuații liniare; această comandă aparține pachetului
linalg, deci înainte de utilizare trebuie încărcat acest pachet sau comanda trebuie
apelată sub forma linalg[linsolve]. Parametrii acestei comenzi sunt: matricea
sistemului A, vectorul termenilor liberi b, și opțional un al treilea parametru, căruia i
se va atribui rangul matricei. Dacă vectorul termenilor liberi se înlocuiește cu o matrice
B, atunci comanda întoarce matricea X care verifică AX = B.
Exemple
> solve(cos(x)+y=9,x);
> solve({cos(x)+y=9, x+y-Pi/2=9},{x,y});
> solve(x^2-2*x+1,x);
> solve(x^2-2*x+7,x);
> dsolve({diff(y(t),t)+t=0, y(1)=2}, y(t));
> isolve(3*x+4*y=7);
> isolve(3*x+4*y=7,t);
> msolve({2*x+3*y=1,5*x-2*y=4},5);
> rsolve(f(n)=2*f(n-1)+3*f(n-2), f(k));
> with(linalg):
> A:=matrix(3,3,[1,2,3,6-8,7,8,1,0,1]);
( )arccos y 9
{ },y 9 x
2
,1 1
,1 6 I 1 6 I
( )y t t2
2
5
2
{ },x 1 4 _Z1 y 1 3 _Z1
{ },x 1 4 t y 1 3 t
{ },x 1 y 3
3
4( )f 0
1
4( )f 1 ( )-1 k
1
4( )f 0
1
4( )f 1 3k
-
Metode Numerice - Lucrarea de laborator 5
3
> b:=vector(3,[-7,2,3]);
> linsolve(A,b);
> linsolve(A,b,'r');
> r;
> B:=matrix(3,2,[-7,3,2,4,3,5]);
> linsolve(A,B);
> C:=matrix(2,2,[1,2,-5,-10]);
> d:=vector(2,[7,-35]);
> linsolve(C,d);
> e:=vector(2,[7,-34]);
> linsolve(C,e)
Ultimul sistem (sistemul Cx = b1) este incompatibil. În această situație MAPLE
întoarce constanta NULL. Această explică lipsa răspunsului la ultima comandă.
:= A
1 2 3
-2 7 8
1 0 1
:= b [ ], ,-7 2 3
, ,
-34
3
-58
3
43
3
, ,
-34
3
-58
3
43
3
3
:= B
-7 3
2 4
3 5
-34
3-2
-58
3-8
43
37
:= C
1 2
-5 -10
:= d [ ],7 -35
[ ],7 2 _t1
_t1
:= e [ ],7 -34
-
Mădălina Roxana Buneci Metode Numerice –Laborator
4
Rezultatele pe care le vom obține în secțiunile următoare pot fi comparate cu
rezultatele comenzii fsolve. Comanda fsolve admite un set de opțiuni. Specificarea lor
nu este obligatorie. Astfel forma generală a comenzii este:
> fsolve(ecuatii, necunoscute (opțional), alți parametrii opționali)
Opțiunile (parametrii opționali) sunt:
complex: caută rădăcinile ecuației în mulțimea numerelor complexe (în virgulă
mobilă).
fulldigits: previne micșorarea numărului de cifre semnificative pe măsură ce se
execută calculele; precizia evaluării crește dar viteză de execuție scade.
optiuni
maxsols = n: calculează cele mai mici n rădăcini; opțiunea are sens doar pentru
ecuațiile polinomiale, pentru că numai pentru acestea se calculează mai multe rădăcini.
interval: caută rădăcinile în intervalul dat; intervalul se specifică sub forma
a..b, sau x = a..b, sau {x = a..b, y = c..d, …}; intervalul se consideră închis. Dacă se
folosește împreună cu opțiunea complex, mulțimea de numere complexe în care se
caută soluția se specifică sub forma unei zone rectangulare din plan dată prin colțul
stânga-jos și colțul dreapta-sus: x = a +I*b..c +d* I (unde a≤c și b≤d).
starting_values: specifică valori inițială pentru o variabilă, sau pentru lista
variabilelor; specificarea se face, de exemplu, sub forma x = a sau {x = a, y = b}, sau
[a, b]. Dacă se utilizează această opțiune nu se mai pot specifica explicit variabilele.
avoid specifică valorile excluse la aplicarea algoritmului de rezolvare
(aproximare). Specificarea se face, de exemplu, sub forma avoid= {x = a, x=b} sau
avoid={ {x = a, y = b},{x=c,y=d} }sau avoid={[a, b]}.
Exemple
> fsolve(tan(sin(x))=1,x,fulldigits);
> fsolve(tan(sin(x))=1,x);
> fsolve(x^4-1,x);
> fsolve(x^4-1,x,0..4);
0.9033391108
0.9033391108
,-1.000000000 1.000000000
-
Metode Numerice - Lucrarea de laborator 5
5
> fsolve(x^4-1,x,complex);
> fsolve(x^4-1,x,complex, maxsols=2);
> Digits:=20;
> fsolve(x^2=2,x=0..infinity);
> fsolve(x^2=2,x=0..infinity,fulldigits);
III.2. Metoda bisecției (metoda înjumătățirii intervalului)
Fie f : [a,b] R, o funcție continuă cu proprietatea că
f(a)f(b) < 0.
Atunci există cel puțin o rădăcină x* (a,b) a ecuației f(x)=0. Pentru găsirea rădăcinii
se micșorează la fiecare pas intervalul în care funcția își schimbă semnul. Metoda
bisecției presupune înjumătățirea la fiecare pas a acestui interval. Astfel
se determină mijlocul c = 2
ba al intervalului (a,b).
dacă f(c)f(a)
-
Mădălina Roxana Buneci Metode Numerice –Laborator
6
Algoritm
Date de intrare:
f continuă, a,b cu f(a)f(b)
-
Metode Numerice - Lucrarea de laborator 5
7
x8 –3x+3 = 0.
Reprezentăm grafic funcția
x->x8 –3x+3
pentru a localiza rădăcinile.
> plot(x^8-3*x-3,x,color=black);
> plot(x^8-3*x-3,x=-5..5,color=black);
> plot(x^8-3*x-3,x=-2..2,color=black);
> plot(x^8-3*x-3,x=-1.5..1.5,color=black);
-
Mădălina Roxana Buneci Metode Numerice –Laborator
8
Se observă că ecuația are două rădăcini reale. Una în intervalul (-1.5, 0) și alta în
intervalul (1,1.5).
> f:=x-> x^8-3*x-3:
> bisectie(f,1,3/2,10^(-3));
> evalf(bisectie(f,1,3/2,10^(-3)));
> bisectie(f,1,1.5,10^(-3));
> bisectie(f,-1.5,0,10^(-3));
> fsolve(f(x),x);
> evalf(bisectie(f,1,3/2,10^(-10)));
> evalf(bisectie(f,1,-3/2,10^(-10)));
III.3. Metoda coardei
Fie f : [a,b] R, o funcție continuă cu proprietatea că
f(a)f(b) < 0.
Rădăcina a ecuației f(x)=0 se caută ca și în cazul metodei bisecției prin micșorarea la
fiecare pas a intervalului în care funcția își schimbă semnul. Metoda coardei presupune:
5207
4096
1.271240234
1.271240234
-0.8801879883
,-0.8800582880 1.271068437
1.271068437
-0.8800582880
-
Metode Numerice - Lucrarea de laborator 5
9
se determină punctul c în care coarda AB intersectează axa Ox, unde A(a,f(a))
și B(b,f(b)).
dacă f(c)f(a)
-
Mădălina Roxana Buneci Metode Numerice –Laborator
10
Dacă m1 > 0, M1 > 0 sunt astfel încât m1 |f’(x)| M1, atunci unica soluției, x*, a
ecuației satisface inegalitățile:
|x*-xn|
1
n
m
xf
|x*-xn| 1nn1
11 xxm
mM
Semnificație geometrică. Fie f : [a, b] R o aplicație de două ori derivabilă
cu f’(x)0, f”(x) 0 oricare ar end if x[a, b], f(a)f(b)
-
Metode Numerice - Lucrarea de laborator 5
11
Cazul 2. f(b)f”(b) 0 (sau echivalent f(a)f”(a) 0) : pentru orice n 1, xn
reprezintă abscisa punctului de intersecție a axei Ox cu coarda ACn-1 unde A(a,f(a)) și
Cn-1(xn-1, f(xn-1)). Așadar avem subcazurile
2.1. f > 0 (f strict descrescătoare)
2.2. f < 0 (f strict crescătoare)
x0 = a x1 x2 b
1.2. f” 0, f(a) 0
a x2 x1 x0 = b
2.1. f” > 0, f(b) 0
-
Mădălina Roxana Buneci Metode Numerice –Laborator
12
Algoritm
Date de intrare:
f de două ori derivabilă pe [a,b], cu f (x) 0, f”(x) 0 pentru orice x și
f(a)f(b)
-
Metode Numerice - Lucrarea de laborator 5
13
dacă f(a) f”(a) 0 atunci
x1: = b ; x2: = a;
cât timp | x1-x2 | execută
x1: = x2;
x2 : = x1 -
f x1
f x1 f b(x1-b),
sfârșit cât timp,
altfel
x1: = a ; x2: = b;
cât timp | x1-x2 | execută
x1: = x2;
x2 : = x1 -
f x1
f x1 f a(x1-a);
sfârșit cât timp,
sfârșit dacă,
x2 reprezintă xN unde N este cel mai mic număr natural cu proprietatea că
xN –xN-1 .
Procedură MAPLE
> mcoarda:= proc(f,a,b,epsilon)
local x1, x2;
if evalf(f(a)*(D@@2)(f)(a))=epsilon do
x1:=x2;
x2:= x1-f(x1)*(x1-b)/(f(x1)-f(b))
end do;
else
x1:=a; x2:=b;
while evalf(abs(x1-x2))>=epsilon do
-
Mădălina Roxana Buneci Metode Numerice –Laborator
14
x1:=x2;
x2:= x1-f(x1)*(x1-a)/(f(x1)-f(a))
end do;
end if;
return x2
end proc;
Aplicăm această procedură pentru determinarea rădăcinilor reale ale ecuației:
x8 –3x+3 = 0.
În secțiunea precedentă rădăcinile au fost localizate în intervalele (+1.5,0) și (1,1.5).
> f:=(x-> x^8-3*x-3);
> mcoarda(f,1,1.5,10^(-3));
> mcoarda(f,-1.5,0,10^(-3));
> fsolve(f(x),x);
> mcoarda(f,1,1.5,10^(-10));
> mcoarda(f,-1.5,0,10^(-10));
III.4. Metoda tangentei
Metoda tangentei este utilizată pentru determinarea unei rădăcini a ecuației f(x)
= 0. Presupunem că f este derivabilă și că derivata nu se anulează. Rădăcina ecuației
este determinată ca limita unui șir. Se pleacă de la un punct x0 dat. Presupunând că s-a
construit termenul xn-1, termenul xn se determină ca fiind abscisa intersecției dintre
tangenta la graficul funcției în xn-1 și axa Ox.
:= f x x8 3 x 3
1.270281421
-0.8741520730
,-0.8800582880 1.271068437
1.271068437
-0.8800582876
-
Metode Numerice - Lucrarea de laborator 5
15
Ecuația tangentei în xn-1 este:
y – f(xn-1) = f (xn-1)(x – xn-1)
Deci intersecția cu axa Ox se află rezolvând sistemul
n-1 n-1 n-1y- f(x ) f (x ) x-x
y 0
În consecință ,
xn = xn-1 -
n 1
n 1
f x
f x
.
Convergența șirului este determinată de termenul inițial x0. Următoarea
teoremă stabilește condiții suficiente pentru convergența metodei tangentei.
Teoremă (Metoda tangentei). Fie f : [a, b] R o aplicație de două ori
derivabilă cu f’(x)0, f”(x) 0 oricare ar end if x[a, b] și f(a)f(b) 0. În plus, oricare ar end if n 1 au
loc următoarele inegalități:
|x* - xn|
1
n
m
xf
xn xn-1
-
Mădălina Roxana Buneci Metode Numerice –Laborator
16
|x* - xn| 21nn1
2 xxm2
M
unde m1 = x [a,b]
f xinf
și M2 = x [a,b]
f xsup
.
Rata convergenței este pătratică.
Semnificație geometrică. Fie f : [a, b] R o aplicație de două ori derivabilă
cu f’(x)0, f”(x) 0 oricare ar end if x[a, b], f(a)f(b) 0. Atunci șirul (xn)n converge la x*
unica soluție a ecuației f(x) = 0. Pentru orice n 1, xn reprezintă abscisa punctului de
intersecție a axei Ox cu tangenta la graficul lui f în punctul de coordonate (xn-1, f(xn-1)).
Deoarece f’ și f” nu se anulează pe [a, b], rezultă că sunt fie strict pozitive fie strict
negative. Așadar avem
Cazul 1. f” > 0 (f strict convexă)
1.1. f’ > 0 (f strict crescătoare)
1.2. f’ < 0 (f strict descrescătoare)
x2 x1 x0
1.1. f’ > 0, f” > 0
-
Metode Numerice - Lucrarea de laborator 5
17
Cazul 2. f” < 0 (f strict concavă)
2.1.f’ > 0 (f strict crescătoare)
2.2. f’ < 0 (f strict descrescătoare)
x0 x1 x2
1.2. f’ 0, f” > 0
x0 x1 x2
2.1. f’ 0, f” 0
x2 x1 x0
2.2. f’ 0, f” 0
-
Mădălina Roxana Buneci Metode Numerice –Laborator
18
Deci pentru aplicarea metodei tangentei în rezolvarea ecuației f(x) = 0 trebuie stabilite
intervalele de monotonie și intervalele de convexitate/concavitate pentru funcția f.
Dacă a și b sunt capetele unui astfel de interval și dacă f(a)f(b)0. Șirul construit rin metoda
tangentei, având termenul inițial x0 converge la unica rădăcină a ecuației f(x) = 0,
situată în intervalul [a, b].
Algoritm
Date de intrare:
f - în condițiile 1.1,1.2,2.1 sau 2.2
x0 - f(x0)f”(x0)>0
0 (precizia –determină condiția de oprire a iterațiilor)
Date de ieșire: xN cu proprietatea că N este cel mai mic număr natural pentru care
xN – xN-1 |2 < .
unde (xn)n este șirul corespunzător metodei tangentei (xN este considerat o aproximație
satisfăcătoare a unicei soluții a ecuației f(x)=0)
x1 := x0;
x2 : = x1 - 1xf
1xf
;
cât timp | x2 – x1 |2 execută
x1 := x2;
x2 : = x1 - 1xf
1xf
;
sfârșit cât timp
Prezentăm în continuare o variantă a acestui algoritm pentru cazul în care f nu verifică
neapărat condițiile suficiente de convergență. Introducem ca dată suplimentară de
intrare numărul maxim de termeni din șir ce urmează a end if calculați (Nmax).
Condiția de oprire se transformă
| xn - xn-1 |2 < sau n > Nmax
x1 := x0;
-
Metode Numerice - Lucrarea de laborator 5
19
x2 : = x1 - 1xf
1xf
;
n : = 1;
cât timp (| x2 – x1 |2 ) și (n Nmax) execută
x1 := x2;
x2 : = x1 - 1xf
1xf
;
n : = n + 1;
sfârșit cât timp
Trebuie verificat la ieșirea din ciclu dacă f(x1) 0. Dacă problema este bine
condiționată, aceasta condiție va asigura acuratețea aproximației.
Proceduri MAPLE
> mtangenta:=proc(f,x0,epsilon)
local x1,x2,df;
df:=D(f);x1:=x0;x2:=x1-f(x1)/df(x1);
while evalf((x2-x1)^2)>=epsilon do
x1:=x2;x2:=x1-f(x1)/df(x1)
end do;
return x2
end proc;
> mtangentaN:=proc(f,x0,epsilon,Nmax)
local x1,x2,n,df;
df:=D(f);x1:=x0;x2:=x1-f(x1)/df(x1);n:=1;
while (evalf((x2-x1)^2)>=epsilon)and (n
-
Mădălina Roxana Buneci Metode Numerice –Laborator
20
Exemple de utilizare a procedurilor
> with(plots):
> plot(exp(x)+2*x+1,x,color=black);
> plot(exp(x)+2*x+1,x=-2..2,color=black);
> f1:=(x->exp(x)+2*x+1);
> mtangenta(f1,0.1,10^(-5));
> fsolve(f1(x),x);
> plot(sin(x)+x-1,x,color=black);
> f2:=(x->sin(x)+x-1);
> mtangenta(f2,1.1,10^(-8));
:= f1 x ex 2 x 1
-0.7388349460
-0.7388350311
:= f2 x ( )sin x x 1
-
Metode Numerice - Lucrarea de laborator 5
21
> fsolve(f2(x),x);
> mtangentaN(f2,1.1,10^(-8),10);
> xN:=mtangentaN(f2,-10.1,10^(-8),10);
> f2(xN);
III.5. Metoda punctului fix
Definiție. Fie (X,d) un spațiu metric și fie f: X X. Funcție f se numește
contracție dacă și numai dacă există q (0,1) astfel încât
d(f(x), f(y) q d(x,y)
pentru orice x,y X
Definiție. Fie f: X X. Punctul X se numește punct fix pentru f dacă f()
= .
0.5109734294
0.5109734294
0.5099954153
0.5109733047
0.5109734294
,Numar de termeni calculati 3
0.5109734294
19.33165959
9.366076806
-4881.864603
-2422.713182
14288.93783
5139.034940
2315.137476
-96996.81022
0.1024903391 108
-0.5376704605 108
,Numar de termeni calculati 10
:= xN -0.5376704605 108
-0.5376704787 108
-
Mădălina Roxana Buneci Metode Numerice –Laborator
22
Teoremă (metoda punctului fix). Fie (X,d) un spațiu metric complet și fie
f:XX o contracție. Atunci există și este unic x* punct fix pentru f. Punctul x* este
limita unui șir construit după cum urmează:
x0 X dat
xn = f(xn-1), n 1.
Mai mult, dacă numărul q (0,1) este astfel încât d(f(x), f(y) qd(x,y)
pentru orice x,y X, atunci pentru orice n1
1. d(x*, xn) q
1 qd(xn, xn-1)
2. d(x*, xn) nq
1 qd(x1, x0)
3. d(x*, xn) qd(x*, xn-1)
Corolar. Fie (E, ) un spațiu Banach (în particular, E = Rm) și S o
submulțime închisă a lui E. Fie f : S S o funcție cu proprietatea că există un număr
q (0,1) este astfel încât f(x) - f(y) qx -y pentru orice x,y S Atunci există și
este unic x* punct fix pentru f și pentru orice x0E, șirul definit recursiv prin
xn = f(xn-1), n 1
converge la x*.
Mai mult, pentru orice n1
1. x*- xn q
1 qxn - xn-1
2. x*- xn nq
1 qx1- x0
3. x*- xn q x*- xn-1
Observație. În general, rata convergenței pentru metoda punctului fix este
liniară.
Algoritm:
Date de intrare:
f (contracție)
-
Metode Numerice - Lucrarea de laborator 5
23
x0 (termenul inițial al șirului)
(precizia ce determină condiția de oprire: se calculează termenii șirului până
la xN cu proprietatea xN-xN-1 ).
Date de ieșire:
xN (aproximație satisfăcătore (determinată de ) pentru punctul fix).
x1:=x0;
x2:=f(x1);
cât timp 1x2x execută
x1:=x2;
x2:=f(x1);
sfârșit cât timp;
La ieșire x2 este aproximație pentru x*, punctul fix al lui f. Faptul că f este contracție
asigură terminarea programului într-un număr finit de pași. Pentru a evita ciclarea în
situația în care un utilizator ar încerca folosirea algoritmului pentru o funcție care nu
este contracție, se poate stabili de la început un număr maxim de pași ce urmează a end
if executați. De asemenea se poate afișa la fiecare pas diferența dintre termenii
consecutivi curenți. (pentru a observa că șirul nu converge dacă f nu e contracție).
Astfel se poate folosi următoarea variantă a algoritmului:
Date de intrare:
f (contracție)
x0 (termenul inițial al șirului)
(precizia)
Nmax (număr maxim de termeni ai șirului ce urmează a end if calculați)
Condiția de oprire: se calculează termenii șirului până la xN cu proprietatea
xN-xN-1 sau N Nmax).
Date de ieșire:
xN (dacă f este contracție, xN este aproximație satisfăcătore - determinată de
și Nmax - pentru punctul fix al lui f).
x1:=x0;
-
Mădălina Roxana Buneci Metode Numerice –Laborator
24
x2:=f(x1);
n:=1;
cât timp ( 1x2x ) sau (n
-
Metode Numerice - Lucrarea de laborator 5
25
Fie f: [1,2] [1,2], definită prin f(x) = 3 1x . Avem f (x) = 1
3 23
1
(x 1)
și deci
x (1,2)
sup | f (x) |
= 1
3 31
4 1
Ca urmare f este contracție pe intervalul [1,2]. Deci soluția ecuației poate end if aflată
ca limita șirului (xn)n, cu
xn = f(xn-1) = 3 n 1x 1 , x0[1,2].
2) Să se rezolve (în mulțimea numerelor reale) ecuația următoare aplicând
metoda punctului fix
x – cos(x) = 0
Notăm g(x) = x – cos(x). Cum cos(x) [-1, 1], soluțiile ecuației x – cos(x) = 0
se găsesc în intervalul [-1, 1]. Deoarece g (x) = 1 + sin(x) 0 pentru x din intervalul
[-1, 1], rezultă că g este strict crescătore pe [-1, 1] și deci ecuația
x – cos(x) = 0
are cel mult o rădăcină. Cum g(0)g(2
)0, ecuația x – cos(x) = 0 are exact o rădăcină
în intervalul [0, 2
]. Fie f: [0,
2
] [0,
2
], f(x) =
1
2( x +cos(x)). Avem f (x) =
1
2(1
– sin(x)) și deci
x (0, 2)
sup | f (x) |
= 1
2 1
Ca urmare f este contracție pe intervalul [0, 2
]. Deci soluția ecuației poate end if aflată
ca limita șirului (xn)n, cu
xn = f(xn-1) =1
2( xn-1 +cos(xn-1)), x0[0,
2
].
Procedura MAPLE de mai jos implementează metoda punctului fix pentru o
contracție pe un interval închis al lui R.
> punctfix:=proc(f,x0, epsilon)
local x1,x2;
x1:=x0;
x2:=f(x1);
-
Mădălina Roxana Buneci Metode Numerice –Laborator
26
while evalf(abs(x2-x1))>=epsilon do
x1:=x2;
x2:=f(x1)
end do;
return x2;
end proc;
Exemplificăm aplicarea procedurii pentru contracțiile
f1: [1,2] [1,2], definită prin f1(x) = 3 x 1 .
f2: [0, 2
] [0,
2
], f2(x) =
1
2( x +cos(x))
comparative cu aplicarea comenzii fsolve.
> f1:=x->(x+1)^(1/3);
> punctfix(f1,1.5,10^(-5));
> fsolve(f1(x)=x,x);
> f2:=x->1/2*(x+cos(x));
> punctfix(f2,1.,10^(-5));
> fsolve(f2(x)=x,x);
Procedura de mai jos implementează metoda punctului fix pentru o contracție pe Rm.
> punctfixM:=proc(f,x0,epsilon,Nmax)
local m,x1, x2,n, norma,i;
m:=nops(x0);
x1:=x0;
x2:=[seq(f[i](op(x1)), i=1..m)];
n:=1;
norma:=0;
for i from 1 to m do norma:=norma + abs(x2[i]-x1[i]) end do;
:= f1 x ( )x 1( )/1 3
1.324719474
1.324717957
:= f2 x 1
2x
1
2( )cos x
0.7390856959
0.7390851332
-
Metode Numerice - Lucrarea de laborator 5
27
while (evalf(norma) >= epsilon) and (n < Nmax) do
x1:=x2;
x2:=[seq(f[i](op(x1)), i=1..m)];
n:=n+1;
norma:=0;
for i from 1 to m do norma:=norma + abs(x2[i]-x1[i]) end do;
end do;
return x2;
end proc;
III.6. Metoda Newton – cazul m-dimensional
Metoda Newton (varianta m - dimensională, m1) este o generalizare a metodei
tangentei. Este o metodă iterativă de rezolvare a unor sisteme de m ecuații și m
necunoscute:
f(x) = 0,
unde f : G Rm, G Rm.
Convenim să notăm cu x1, x2,…, xn,… un șir de elemente din Rm. Rezervăm
indicii inferiori pentru a desemna componentele unui element x = (x1, x2,…,xm) din Rm.
Dacă G este o mulțime deschisă și f : G Rm este o funcție diferențiabilă pe G,
identificăm diferențiala de ordinul I a lui f în x, df(x), cu jacobianul lui f în x, notat
Jf(x):
Jf(x) = mj,i1
j
i xx
f
În cele ce urmează presupunem că matricea Jf(x) este inversabilă pentru x G.
Metoda Newton constă în aproximarea soluției ecuației considerate cu xn,
unde
xn = xn-1 – Jf(xn-1)-1 f(xn-1) (*)
iar aproximația inițială x0G este suficient de apropiată de soluția sistemului.
Presupunem că f este de clasă C2 și sistemul f(x) = 0 admite o soluție x* G cu
proprietatea că Jf(x*) este matrice inversabilă (sau echivalent det(Jf(x*)) 0). Atunci
-
Mădălina Roxana Buneci Metode Numerice –Laborator
28
există o vecinătate deschisă V G a lui x* astfel încât pentru orice x V să avem Jf(x)
este inversabilă. Considerăm funcția
g: V Rm, definită prin g(x) = x – Jf(x)-1 f(x).
Avem
g(x*) = x* - Jf(x*)-1 f(x*) = x*,
deci x* punct fix pentru g. Cum
Jg (x*) = 0 1,
rezultă că există r 0 astfel încât astfel dacă notăm
S = r*,xB = x Rm, x-x* r
avem S V și gS : S S este contracție. Ca urmare șirul definit prin xn= g(xn-1), n1
converge la x* pentru orice termen inițial x0S. Din definiția lui g rezultă că
xn = xn-1 – Jf(xn-1)-1 f(xn-1), n 1,
adică șirul dat de relația (*) (corespunzător metodei lui Newton pentru f(x) = 0 cu
termenul inițial x0 S).
Observaţie Amplificând relația (*) cu Jf(xn-1) rezultă
Jf(xn-1)(xn - xn-1) = –f(xn-1) (**)
sau echivalent
n 1m i n n 1 n 1
j j ijj 1
f xx x f x
x
, i = 1,2,…, m
Dacă se folosește relația (*) pentru determinarea lui xn este necesar să se calculeze
inversa matricei Jf(xn-1). Dacă se folosește relația (**), este necesar să se rezolve un
sistem liniar cu m ecuații, și necunoscutele n 1 n n 1k k kx x x , k = 1, …, m.
Rezolvarea acestui sistem necesită un număr mai mic de operații decât inversarea
matricei Jf(xn-1). Folosim relația (**) se înlocuiește rezolvarea sistemului neliniar prin
rezolvarea succesivă a unor sisteme liniare.
Metoda Newton este o metodă frecvent folosită deoarece este foarte rapid
convergentă (rata convergenței este pătratică). Dar ca și în cazul metodei tangentei
convergența metodei depinde de alegerea aproximației inițiale. Aproximația inițială
-
Metode Numerice - Lucrarea de laborator 5
29
trebuie luată cât mai aproape de soluția problemei, eventual utilizând o altă metodă de
găsire a soluției.
Parametrii procedurii mnewton (de mai jos) sunt
funcția f (se presupune că se rezolvă sistemul f(x) = 0)
x0 = termenul inițial din șirul definit de (**)
epsilon = precizia
Nmax = numărul maxim de termeni din șir ce vor end if calculați
Se calculează n termeni, cu n verificând
(2
n n 1x x < epsilon) sau (nNmax).
Comanda
>subs(expr1,expr2);
substituie subexpresia expr1 în expresia expr2. Comanda
>jacobian(f(x,y,...), [x, y, ...]);
calculează jacobianul lui f. Este o comandă ce aparține pachetului linalg. Comanda
>norm(a, t);
calculează norma t, unde t =1, 2, infinity a vectorului (sau matricei) a. Este de asemenea
o comanda ce aparține pachetului linalg.
Procedură MAPLE
> mnewton := proc(f, x0, epsilon, Nmax)
local m,x1, x2, dx, b, fx, fx1, n, i, j, ex, r;
uses linalg;
m:=vectdim(x0);x1 := vector(m);x2:=vector(m);
for i to m do x1[i] := x0[i] end do;
dx := vector(m);
b := vector(m);
fx := jacobian(f(seq(x[i],i=1..m)), [seq(x[i],i=1..m)]);
fx1 := matrix(m, m);
ex := seq(x[i] = x1[i], i = 1 .. m);
for i to m do
for j to m do fx1[i, j] := evalf(subs(ex, fx[i, j])) end do
end do;
b := map(-evalf,f(seq(x1[i],i=1..m)));
dx := linsolve(fx1, b, 'r');
if r m then print(`Metoda nu se aplica`); RETURN(NULL) end if;
-
Mădălina Roxana Buneci Metode Numerice –Laborator
30
for i to m do x2[i] := x1[i] + dx[i] end do;
n := 1;
print(x2);
while epsilon f:=(x,y)->[x^2-y,x^3-5*y];
> mnewton(f,vector([10,0.1]),10^(-5),9);
:= f ( ),x y [ ],x2 y x3 5 y
[ ],7.500000001 50.00000002
[ ],6.000000002 33.75000004
[ ],5.250000001 27.00000001
[ ],5.021739128 25.16576084
[ ],5.000186603 25.00140155
[ ],5.000000017 25.00000013
,Numar de pasi 6
,Valoarea functiei [ ],0.00046451 0.0069879
-
Metode Numerice - Lucrarea de laborator 5
31
> fsolve({f(x,y)[1],f(x,y)[2]},{x=10,y=0.1});
> f1:=(x,y)->[x^2+y^2-1,x^3-y];
> mnewton(f1,vector([0.9,0.5]),10^(-5),9);
> fsolve({f1(x,y)[1],f1(x,y)[2]},{x=0.9,y=0.5});
> mnewton(f1,vector([1,1]),10^(-5),9);
> fsolve({f1(x,y)[1],f1(x,y)[2]},{x=1,y=1});
> mnewton(f1,vector([-1,1]),10^(-5),20);
[ ],5.000000017 25.00000013
{ },x 5.000000000 y 25.00000000
:= f1 ( ),x y [ ], x2 y2 1 x3 y
[ ],0.8316784870 0.5629787234
[ ],0.8260617824 0.5636079087
[ ],0.8260313586 0.5636241619
,Numar de pasi 3
,Valoarea functiei [ ],0.000031943 0.0000785347
[ ],0.8260313586 0.5636241619
{ },x 0.8260313577 y 0.5636241622
[ ],0.8750000000 0.6250000000
[ ],0.8290363483 0.5643491124
[ ],0.8260401080 0.5636197732
,Numar de pasi 3
,Valoarea functiei [ ],0.005791188 0.0054486200
[ ],0.8260401080 0.5636197732
{ },x 0.8260313577 y 0.5636241622
[ ],-0.2500000001 1.250000000
[ ],-81.50000065 -15.25000014
[ ],-54.33007595 64.91769907
[ ],-36.21723980 24.89070061
[ ],-24.14389039 3.68381152
[ ],-16.10992954 -24.48752262
[ ],-10.75342471 -10.48890844
[ ],-7.178878965 -3.444508772
-
Mădălina Roxana Buneci Metode Numerice –Laborator
32
> fsolve({f1(x,y)[1],f1(x,y)[2]},{x=-1,y=1});
> f2:=(x,y,z)->[x+x^2-2*y*z-0.1,y-y^2+3*x*z+0.2, z+z^2+2*x*y-0.3];
> mnewton(f2,vector([0,0,0]),10^(-5),10);
> fsolve({f2(x,y,z)[1],f2(x,y,z)[2], f2(x,y,z)[3]},{x=0,y=0,z=0});
[ ],-4.781922755 0.617903126
[ ],-3.270919075 -5.691810662
[ ],-2.260819460 -2.574374564
[ ],-1.578178423 -1.088178100
[ ],-1.121843867 -0.5209805009
[ ],-0.8857962879 -0.5206567978
[ ],-0.8296024170 -0.5627517822
[ ],-0.8260437059 -0.5636179649
[ ],-0.8260313579 -0.5636241623
,Numar de pasi 17
,Valoarea functiei [ ],0.000013414 -0.0000314743
[ ],-0.8260313579 -0.5636241623
{ },x -0.8260313577 y -0.5636241622
:= f2 ( ), ,x y z [ ], , x x2 2 y z 0.1 y y2 3 x z 0.2 z z2 2 x y 0.3
[ ], ,0.1 -0.2 0.3
[ ], ,0.02245322250 -0.1743243244 0.2461538462
[ ], ,0.01287849239 -0.1778109522 0.2447473526
[ ], ,0.01282415092 -0.1778006638 0.2446880471
,Numar de pasi 4
,Valoarea functiei [ ], ,0.0000818676 0.0000282439 0.0000687452
[ ], ,0.01282415092 -0.1778006638 0.2446880471
{ }, ,x 0.01282414583 y -0.1778006680 z 0.2446880443