Curs 3 MN Ecuatii neliniare.pdf

CURSUL NR.3 Metode numerice Prof. Dumitru Nicoar

1

Rezolvarea numeric a ecuaiilor algebrice i transcendente

Exist numeroase situaii n care avem de rezolvat ecuaii polinomiale sau

nepolinomiale (transcendente) cu o singur variabil de forma

f ( x ) = 0 (1)

ale cror soluii nu se pot determina prin metodele algebrice cunoscute. Pentru

rezolvarea acestor ecuaii este necesar mai nti s se identifice, printr-o anumit

metod, intervalele n care se afl exact o rdcin a ecuaiei.

1. Metoda grafic

O metoda simpl pentru a obine o estimare (aproximare) a

radaciniilor unei ecuaii f(x)=0 este sa facem graficul

funciei si sa observam unde i de cte ori graficul

intersecteaza axa x.

Fcnd graficul funciei putem avea urmatoarele situaii:

a) acelai semn - nu sunt rdcini

b) semne diferite o singr

c) acelai semn - dou rdcini

d) semne diferite- trei radacini

a b


2

Metoda grafic utila pentru a afla daca exista si care este numarul rdcinilor

Aplicaii metoda grafic: radcinile ecuaiei f(x)=0

Exemplu 1.

0sin xx sau xx sin (2)

Facem graficele pentru:

y1=sin(x) i y2=x; (3)

!!!! Numarul radcinilor coincide cu numarul punctelor de interseci ale graficelor

>>x=0:0.01:2*pi;

>> y1=sin(x);

>>plot(x,y1)

>> hold on

>> y2=x;

>> plot(x,y2)

>> grid

>> gtext('x')

>> gtext('sin x')

Figura 1

Observam ca exista un singur punct comun in x=0, care este singura radacina.

0 1 2 3 4 5 6 7-1

0

1

2

3

4

5

6

7

sin x

x


3

0 2 4 6 8 10 12 14-10

-8

-6

-4

-2

0

2

4

6

8

10

x

tan x

tan x

tan x

pi/2 3*pi/2

Exemplu 2.

0tan xx sau xx tan (4)

>>x=0:0.01:4*pi;

>> y1=tan(x);

>> plot(x,y1)

>> y2=x;

>> hold on

>> plot(x,y2)

>> grid

>> gtext('x')

>> gtext('tan x')

Figura 2

Ecuatia are o infinitate de solutii: (x = 0,4.493,7.725, . . .).

Cum determinm soluiile?

Teorema 1.1 (Teorema valorii intermediare) Dac funcia f(x) este continu n

intervalul [a, b] atunci pentru orice y dintre f(a) i f(b) exista ],[ bac astel nct

ycf )( .

Teorema este util n condiiile urmtoare: funcia f(x) este continu n intervalul

[a, b] , strict monoton, adic derivata f '(x) are semn constant pe intervalul [a,b]

i valorile funciei n a i b, f(a) i f(b) au semne diferite adica f(a) * f(b) < 0 ,

atunci exist o singur radcin a ecuaiei 0)( xf n intervalul [a, b].

Rbaf ],[: continu i strict monoton;

;0)()( bfaf

xf are o singur rdcin pe [a,b].


4

Metodele cele mai utilizate n calculul numeric al rdcinilor unei ecuaii care

satisface ipotezele de mai sus sunt:

1. Metoda njumtirii intervalului - metoda biseciei;

2. Metoda secantei;

3. Metoda Newton (Raphson) (metoda tangentelor).

1. Metoda njumtirii intervalului (biseciei)

Metoda biseciei sau metoda njumtirii intervalului este una din primele metode

dezvoltate pentru determinarea rdacinilor ecuaiei neliniare f(x) = 0 .

Algoritmul metodei:

pasul1: iniializm numarul iteraiei k=1;

pasul 2: fie 2

bac

dac 0)( cf sau 02

ab

stop ;

pasul3: dac

0)()( cfaf ca atunci rdcina r

x este n intercalul [c,b]

altfel

0)()( cfaf cb atunci rdcina r

x este n intercalul [a,c]

(cazul din figura 3)

i napoi la pasul 2.

Figura 3


5

Mrimea intervalului [an, bn ] dup n pai este

n

nn

nn

ababbaI

22],[

(5)

Mrimea intervalului la fiecare iteraie corespunde cu maximul erorii pentru r

x ,

astfel

nn

abE

2

(6)

Dac dorim ca eroarea s fie mai mic dect

n

ab

2

abn 2 2ln

ab

n

2ln443.1

abn

(7)

Observaie. Trebuie s executm n iteraii pentru a avea ca ordin de

convergen.

Algoritmul este implementat in m-fila bisection_ro.m

function [x,k] = bisection_ro(f,a,b,tol,M,display)

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % % % BISECTION.M Routina pentru implementarea metodei bisectiei % % % %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % % Inputs: f = functia se da fie ca annimus function fie dintr-o % m-fila tip function de forma: function y = f(x) % a,b = capetele intervalului(unde b > a) % tol = criteriu de convergenta (tol > 100*eps) % M = maximum numarul maxim de iteratii (M >= 2) % dispaly = rezultatele intermediare sunt afisate daca > 0

. % % Outputs: x = radacina estimata a ecuatiei f(x) = 0 % k = numarul de iteratii efectuate % % Observatie. a si b trebuie ales astfel ca f(a)*f(b) < 0 pentru a fi


6

% siguri ca in intervalul [a,b] este o radacina % % %%%%% Fila adaptata dupa J. R. White, UMass-Lowell (July 2003)%%%%%%% % function [x,k] = bisection_ro(f,a,b,tol,M,display) % % cerceteaza marginile, toleranta, # iteratiile, si display valorile

% intermediare

% calculeaza valorile functiei cu m-fila feval.m (help feval) fa = feval(f,a); fb = feval(f,b); if fa*fb >= 0 fprintf ('\n !Atentie ! In bisection.m trebuie f(a)f(b) < 0.\n\n'); return end if tol < 100*eps, tol = 100*eps; end if M < 2, M = 2; end if display > 0 fprintf('\n\n Valorile intermediare pentru metoda bisectiei: \n\n') fprintf(' k a f(a) b f(b) xr ...

f(xr) semn(fa*fr)\n') end

% determina radacinile k = 1; err = tol + 1; while (err > tol) & (k 0 fprintf(' %3i %10.6f %10.6f %10.6f %10.6f %10.6f %10.6f %6i \n',

... k,a,fa,b,fb,xr,fr,sign(fa*fr)); end if fa*fr < 0 b = xr; fb = fr; else a = xr; fa = fr; end err = abs(fr); k = k + 1; end x = xr; k = k-1; % save outputs if k == M fprintf('\n\n Aten?ie: Numrul maxim de itera?ii a fost atins!!! \n') end % % end of function

Observaie. eps - precizia relativ n virgul mobil (help eps)

Exemplul 1

23)(3 xxxf

Functia este dat ca functie anonimus @(x) x^3-3*x-2


7

>> bisection(@(x)x^3-3*x-2, 0, 4, 200*eps, 20, 1)

Valorile intermediare pentru metoda bisectiei:

k a f(a) b f(b) xr f(xr) sign(fa*fr)

1 0.0000 -2.0000 4.0000 50.0000 2.0000 0.0000 0

ans =2

Exemplul 2

7)(2 xxf

>> bisection_ro(@(x) x.^2-7, 1, 4, 100*eps, 10, 1)


k a f(a) b f(b) xr f(xr) semn(fa*fr)

1 1.000000 -6.000000 4.000000 9.000000 2.500000 -0.750000 1

2 2.500000 -0.750000 4.000000 9.000000 3.250000 3.562500 -1

3 2.500000 -0.750000 3.250000 3.562500 2.875000 1.265625 -1

4 2.500000 -0.750000 2.875000 1.265625 2.687500 0.222656 -1

5 2.500000 -0.750000 2.687500 0.222656 2.593750 -0.272461 1

6 2.593750 -0.272461 2.687500 0.222656 2.640625 -0.027100 1

7 2.640625 -0.027100 2.687500 0.222656 2.664063 0.097229 -1

8 2.640625 -0.027100 2.664063 0.097229 2.652344 0.034927 -1

9 2.640625 -0.027100 2.652344 0.034927 2.646484 0.003880 -1

10 2.640625 -0.027100 2.646484 0.003880 2.643555 -0.011619 1

Atentie: Numrul maxim de iteratii a fost atins !!!

ans = 2.6436


8

Exemplu 3.

xxxf tan)(

Varianta in care se da functia ca anonimus function @(x) tan(x)-x

>> bisection(@(x)tan(x)-x,2,4.6,200*eps,50,1)


k a f(a) b f(b) xr f(xr) sign(fa*fr)

1 2.000000 -4.185040 4.600000 4.260175 3.300000 -3.140254 1

2 3.300000 -3.140254 4.600000 4.260175 3.950000 -2.902889 1

3 3.950000 -2.902889 4.600000 4.260175 4.275000 -2.136396 1

4 4.275000 -2.136396 4.600000 4.260175 4.437500 -0.891762 1

5 4.437500 -0.891762 4.600000 4.260175 4.518750 0.580791 -1

6 4.437500 -0.891762 4.518750 0.580791 4.478125 -0.287812 1

7 4.478125 -0.287812 4.518750 0.580791 4.498437 0.103984 -1

8 4.478125 -0.287812 4.498437 0.103984 4.488281 -0.101095 1

9 4.488281 -0.101095 4.498437 0.103984 4.493359 -0.001011 1

10 4.493359 -0.001011 4.498437 0.103984 4.495898 0.050851 -1

11 4.493359 -0.001011 4.495898 0.050851 4.494629 0.024764 -1

12 4.493359 -0.001011 4.494629 0.024764 4.493994 0.011838 -1

13 4.493359 -0.001011 4.493994 0.011838 4.493677 0.005404 -1

14 4.493359 -0.001011 4.493677 0.005404 4.493518 0.002194 -1

15 4.493359 -0.001011 4.493518 0.002194 4.493439 0.000591 -1

16 4.493359 -0.001011 4.493439 0.000591 4.493399 -0.000210 1

17 4.493399 -0.000210 4.493439 0.000591 4.493419 0.000190 -1

18 4.493399 -0.000210 4.493419 0.000190 4.493409 -0.000010 1

19 4.493409 -0.000010 4.493419 0.000190 4.493414 0.000090 -1

20 4.493409 -0.000010 4.493414 0.000090 4.493411 0.000040 -1

21 4.493409 -0.000010 4.493411 0.000040 4.493410 0.000015 -1

22 4.493409 -0.000010 4.493410 0.000015 4.493410 0.000003 -1

23 4.493409 -0.000010 4.493410 0.000003 4.493409 -0.000004 1

24 4.493409 -0.000004 4.493410 0.000003 4.493409 -0.000001 1

25 4.493409 -0.000001 4.493410 0.000003 4.493410 0.000001 -1

26 4.493409 -0.000001 4.493410 0.000001 4.493409 0.000000 -1

ans = 4.4934


9

Exemplu 3.

xxxf tan)(

Varianta in care se da functia printr-o m-fila function function f = ftg(x) f=tan(x)-x;

>> [xr,k] = bisection_ro('ftg',2,4.6,1e-6,50,1)


k a f(a) b f(b) xr f(xr) semn(fa*fr)

1 2.000000 -4.185040 4.600000 4.260175 3.300000 -3.140254 1

2 3.300000 -3.140254 4.600000 4.260175 3.950000 -2.902889 1

3 3.950000 -2.902889 4.600000 4.260175 4.275000 -2.136396 1

4 4.275000 -2.136396 4.600000 4.260175 4.437500 -0.891762 1

5 4.437500 -0.891762 4.600000 4.260175 4.518750 0.580791 -1

6 4.437500 -0.891762 4.518750 0.580791 4.478125 -0.287812 1

7 4.478125 -0.287812 4.518750 0.580791 4.498437 0.103984 -1

8 4.478125 -0.287812 4.498437 0.103984 4.488281 -0.101095 1

9 4.488281 -0.101095 4.498437 0.103984 4.493359 -0.001011 1

10 4.493359 -0.001011 4.498437 0.103984 4.495898 0.050851 -1

11 4.493359 -0.001011 4.495898 0.050851 4.494629 0.024764 -1

12 4.493359 -0.001011 4.494629 0.024764 4.493994 0.011838 -1

13 4.493359 -0.001011 4.493994 0.011838 4.493677 0.005404 -1

14 4.493359 -0.001011 4.493677 0.005404 4.493518 0.002194 -1

15 4.493359 -0.001011 4.493518 0.002194 4.493439 0.000591 -1

16 4.493359 -0.001011 4.493439 0.000591 4.493399 -0.000210 1

17 4.493399 -0.000210 4.493439 0.000591 4.493419 0.000190 -1

18 4.493399 -0.000210 4.493419 0.000190 4.493409 -0.000010 1

19 4.493409 -0.000010 4.493419 0.000190 4.493414 0.000090 -1

20 4.493409 -0.000010 4.493414 0.000090 4.493411 0.000040 -1

21 4.493409 -0.000010 4.493411 0.000040 4.493410 0.000015 -1

22 4.493409 -0.000010 4.493410 0.000015 4.493410 0.000003 -1

23 4.493409 -0.000010 4.493410 0.000003 4.493409 -0.000004 1

24 4.493409 -0.000004 4.493410 0.000003 4.493409 -0.000001 1

xr = 4.49340943098068 k = 24


10

Observaie. Metoda biseciei eueaz dac:

Funcia are mai mult de o rdcin n intervalul [a,b];

Funcia are rdcin dubl;

Funcia are singulariti (este discontinu ntr-un punct din intervalul in care

se caut rdcina).

2. Metoda secantei

Se tie din cele de mai sus c pentru o funcia f(x) este continu n intervalul [a, b]

i strict monoton, adic derivata f '(x) are semn constant pe intervalul [a,b] i

valorile funciei n a i b, f(a) i f(b) au semne diferite adica f(a) * f(b) < 0 , atunci

exist o singur radcin a ecuaiei 0)( xf n intervalul [a, b]. Aceasta rdcin

se poate aproxima cu abscisa punctului de intersece a secantei (coardei) care trece

prin punctele A(a,f(a)) i B(b,f(b)) cu axa Ox, figura 4.

Figura 4

Ecuaia secantei este:

axab

afbfafy

)()()( (8)

Din (8) se obine abscisa

)()(

)(1

afbf

abafax

(9)

dac

0)()(1 xfaf noul subinterval este [a, x1] (ca n figura 4).

altfel


11

0)()(1 xfaf noul subinterval este [x1, b] i

procedeul se repet dupa relaia de recuren

)()(

)(1

1

1

nn

nn

nnnxfxf

xxxfxx . (10)

Algoritmul este implementat in m-fila secant_ro.m

function [x,k] = secant(f,x0,x1,tol,M,display)

%

% SECANT.M Routina pentru implementarea metodei secantei % % Inputs: f = functia % The function f is of the form: function y = f(x) % x0,x1 = valorile initiale intre care este localizata solutia % tol = criteriul de convergenta (tol > 100*eps) % M = numarul maxim de iteratii (M >= 2) % display = rezultatele intermediare sunt afisate pentru valoare % dysplay > 0 . % % Outputs: x = estimarile pentru radacinile f(x) = 0 % k = numarul iteratiilor % % Fila adaptata dupa J. R. White, UMass-Lowell (July 2003) % function [x,k] = secant(f,x0,x1,tol,M,display) % % check tolerance, # iterations, initial guesses, and display switch if tol < 100*eps, tol = 100*eps; end if M < 2, M = 2; end if abs(x1-x0) < tol fprintf ('The Secant Method requires that x1 be different from

x0.\n'); return end if display > 0 fprintf('\n\n Intermediate edit from the Secant Method: \n\n'); fprintf(' k x0 f0 x1 f1 x2

f2\n'); end % % find root f0 = feval(f,x0); f1 = feval(f,x1); k = 1; err = tol + 1; while (err > tol) & (k 0 fprintf(' %3i %10.6f %10.6f %10.6f %10.6f %10.6f %10.6f \n', ... k,x0,f0,x1,f1,x2,f2); end err = abs(f2); k = k + 1;


12

x0 = x1; f0 = f1; x1 = x2; f1 = f2; end x = x2; k = k-1; % save outputs if k == M fprintf('\n\n Warning: The maximum number of iterations was

reached!!! \n') end % % end of function

Se apeleaz cu comanda

>> secant(f,x0,x1,tol,M,display) Exemplu 2.1

>> secant(@(x) log(x+1)+6*x.^2-3*x-1,1,3,200*eps,20,1)

Intermediate edit from the Secant Method:

k x0 f0 x1 f1 x2 f2

1 1.000000 2.693147 3.000000 45.386294 0.873837 1.588024

2 3.000000 45.386294 0.873837 1.588024 0.796747 1.004573

3 0.873837 1.588024 0.796747 1.004573 0.664016 0.162688

4 0.796747 1.004573 0.664016 0.162688 0.638366 0.023670

5 0.664016 0.162688 0.638366 0.023670 0.633999 0.000762

6 0.638366 0.023670 0.633999 0.000762 0.633854 0.000004

7 0.633999 0.000762 0.633854 0.000004 0.633853 0.000000

8 0.633854 0.000004 0.633853 0.000000 0.633853 0.000000

ans = 0.6339

Exemplu 2.2

>> secant(@(x) tan(x)-x,2,4,200*eps,50,1)

Intermediate edit from the Secant Method:

k x0 f0 x1 f1 x2 f2

1 2.000000 -4.185040 4.000000 -2.842179 8.233020 -10.743705

2 4.000000 -2.842179 8.233020 -10.743705 2.477383 -3.260255


13

3 8.233020 -10.743705 2.477383 -3.260255 -0.030129 -0.000009

4 2.477383 -3.260255 -0.030129 -0.000009 -0.030136 -0.000009

5 -0.030129 -0.000009 -0.030136 -0.000009 -0.020091 -0.000003

6 -0.030136 -0.000009 -0.020091 -0.000003 -0.015862 -0.000001

7 -0.020091 -0.000003 -0.015862 -0.000001 -0.011765 -0.000001

8 -0.015862 -0.000001 -0.011765 -0.000001 -0.008941 -0.000000

9 -0.011765 -0.000001 -0.008941 -0.000000 -0.006732 -0.000000

10 -0.008941 -0.000000 -0.006732 -0.000000 -0.005087 -0.000000

11 -0.006732 -0.000000 -0.005087 -0.000000 -0.003839 -0.000000

12 -0.005087 -0.000000 -0.003839 -0.000000 -0.002898 -0.000000

13 -0.003839 -0.000000 -0.002898 -0.000000 -0.002188 -0.000000

14 -0.002898 -0.000000 -0.002188 -0.000000 -0.001651 -0.000000

15 -0.002188 -0.000000 -0.001651 -0.000000 -0.001247 -0.000000

16 -0.001651 -0.000000 -0.001247 -0.000000 -0.000941 -0.000000

17 -0.001247 -0.000000 -0.000941 -0.000000 -0.000710 -0.000000

18 -0.000941 -0.000000 -0.000710 -0.000000 -0.000536 -0.000000

19 -0.000710 -0.000000 -0.000536 -0.000000 -0.000405 -0.000000

20 -0.000536 -0.000000 -0.000405 -0.000000 -0.000306 -0.000000

21 -0.000405 -0.000000 -0.000306 -0.000000 -0.000231 -0.000000

22 -0.000306 -0.000000 -0.000231 -0.000000 -0.000174 -0.000000

23 -0.000231 -0.000000 -0.000174 -0.000000 -0.000131 -0.000000

24 -0.000174 -0.000000 -0.000131 -0.000000 -0.000099 -0.000000

25 -0.000131 -0.000000 -0.000099 -0.000000 -0.000075 -0.000000

26 -0.000099 -0.000000 -0.000075 -0.000000 -0.000057 -0.000000

27 -0.000075 -0.000000 -0.000057 -0.000000 -0.000043 -0.000000

ans = -4.268204728757503e-005

3. Metoda Newton-Raphson

3.1 Algoritmul metodei

Metoda Newton-Raphson poate fi folosit pentru a rezolva problemele de mai sus,

n care metoda biseciei sau secantei nu se pot aplica. Metoda cere ca s existe

prima derivat a funciei f(x) i s fie continu n vecintatea soluiei.

Strategia metodei Newton-Raphson const n aproximarea curbei f(x) prin tangenta

sa ntr-o estimare oarecare, de exemplu n estimarea xk obnut la pasul la pasul k,

figura 5. Ecuaia tangentei (o dreapt) n punctul ))(,( kk xfx este:

)()()( kkk xxxfxfy (11)


14

La pasul urmtor, alegem estimarea xk+1 punctul n care tangenta intersecteaz axa

ox

)()()(0 1 kkkk xxxfxf (12)

de unde rezult formula iterativ a metodei Newton:

)(

)(1

k

kkk

xf

xfxx

(13)

Figura 5

Algoritmul Newton- Raphson este implementat n rutina MATLAB newton.m

function [x,fx,xx] = newton(f,df,x0,TolX,MaxIter) %newton.m pentru solutia f(x) = 0 folosind metoda Newton. % %input: f = functia care se va defini printr-o M-fila sau inline % df = df(x)/dx (daca nu este data derivata de utilizator) % x0 = valoarea initiala de start al algoritmului % TolX = limita superioara pentru |x(k) - x(k-1)| % MaxIter = maximum pentru numarul de iteratii % %output: x = valoarea (punctul) obtinuta de algoritm % fx = f(x(last)), xx = istoria lui x % h = 1e-4; h2 = 2*h; TolFun=eps; if nargin == 4 & isnumeric(df), MaxIter = TolX; TolX = x0; x0 = df; end xx(1) = x0; fx = feval(f,x0); for k = 1: MaxIter if ~isnumeric(df), dfdx = feval(df,xx(k)); %derivative function


15

else dfdx = (feval(f,xx(k) + h)-feval(f,xx(k) - h))/h2; %numerical drv end dx = -fx/dfdx; xx(k+1) = xx(k)+dx; %Eq.(4.4.2) fx = feval(f,xx(k + 1)); if abs(fx)>f1 = inline('tan(pi - x)-x','x');

>>x0 = 1.8; TolX = 1e-5; MaxIter = 50; % cu valoarea de start 1.8

>>[x,err,xx] = newton(f1,x0,1e-5,50)

b) derivata este dat de utilizator ( mai jos este dat prin comanda inline

>>df1 = inline('-(sec(pi-x)).^2-1', 'x'); % prima derivat

>>[x,err,xx1] = newton(f1,df1,1.8,1e-5,50)

Caz a)

>> f1 = inline('tan(pi - x)-x','x');

>> x0 = 1.8; TolX = 1e-5; MaxIter = 50;

>> [x,err,xx] = newton(f42,1.8,1e-5,50)

x = 2.0288

err =1.1684e-011

xx =

1.8000 1.9220 2.0075 2.0280 2.0288 2.0288

3.2 Analiza erorii pentru metoda Newton

Pentru a face o analiz a erorilor vom considera dezvoltarea taylor de ordinul doi a

funciei f(x) n punctul x = xk:


16

2)(2

)()()()()( k

kkkk xx

xfxxxfxfxf

(14)

Vom substitui x cu soluia 0x , 0)( 0 xf i vom obine:

200 )(

2

)()()()( k

kkkk xx

xfxxxfxf

(15)

Itroducem (15) n relaia (13) si obinem:

2001 )(

)(2

)()( k

k

kkkk xx

xf

xfxxxx

(16)

Definim eroarea estimrii xk prin 0xxe kk i obinem

kkkkkk

k

kk eeAeAe

xf

xfe

221

)(2

)( (17)

Din relaia (17) rezult c dac marimea eroarii iniiale 0e este mic, astfel ca ,

10 Ae , atunci i ordinul de mrime al erorilor succesive va fii mic, asta dac Ak

nu devine prea mare.

Spunem c metoda Newton este ptratic convergent deoarece, aa cum se vede

din relaia (17), marimea eroarii estimate la un anumit pas este proporional cu

ptratul erorii estimate la pasul anterior.

Observaie asupra metodei Newton:

Converge foarte repede;

Metoda poate determina i rdcinile complexe.

Metoda eueaz dac:

Punctul inial este un punct de min/max local (tangenta nu mai intersecteaz

axa Ox);

Punctul inial este ales greit datorit formei graficului;

Curs 3 MN Ecuatii neliniare.pdf

Documents

Transcript of Curs 3 MN Ecuatii neliniare.pdf