6a lezione - laboratorio

a.a 2004-2005

Corso di Laurea ING. MECCANICA

La La scatolascatola deidei cocoloriloriAvevo una scatola di coloriAvevo una scatola di coloriOgni colore splendeva di gioiaOgni colore splendeva di gioiaAvevo una scatola di coloriAvevo una scatola di coloriCaldi e delicati e vivi.Caldi e delicati e vivi.Non avevo il rosso per le ferite e per il sangue,Non avevo il rosso per le ferite e per il sangue,Non avevo il nero per l’orfano,Non avevo il nero per l’orfano,Non avevo il bianco per il viso dei morti,Non avevo il bianco per il viso dei morti,Non avevo il giallo per le sabbie infuocate,Non avevo il giallo per le sabbie infuocate,Avevo l’arancione per la gioia e per la vita,Avevo l’arancione per la gioia e per la vita,Avevo il verde per germogli e fioriture,Avevo il verde per germogli e fioriture,Avevo il blu per limpidi cieli azzurri,Avevo il blu per limpidi cieli azzurri,Avevo il rosa per i sogni e per il riposo.Avevo il rosa per i sogni e per il riposo.Mi sedetti Mi sedetti ee dipinsi dipinsiLa PACE La PACE Tali ShurekTali Shurek

Esercizio 1

Si verifica facilmente che x = 0 è soluzione.b) Senza preoccuparsi delle ipotesi di convergenza, si approssimi tale soluzione, applicando il metodo di Newton e quello di bisezione per 15 iterazioni.c) Fare le dovute considerazioni.

Data l’equazione:

sin 0.5sin 2 0, 1, 1f x x x x

Punto b: metodo di Newton x0=-1;nmax=15; toll=1e-12;fun='sin(x)-0.5*sin(2*x)';dfun='cos(x)-cos(2*x)';[xvect,xdiff,fx,it,p,c]=newton(x0,nmax,toll,fun,dfun); Superato il numero massimo di iterazioni Numero di Iterazioni : 15 Radice calcolata : -1.9309205276636407e-003 Ordine stimato : 0.9999897626825431 Fattore di riduzione : 0.6666198373211456

Tabella riassuntiva - Newton iter=0:it;fprintf('%2d %23.15e %15.3e %15.3e\n',[iter' xvect xdiff fx]')

0 -1.000000000000000e+000 0.000e+000 3.868e-001 1 -5.955642027988763e-001 4.044e-001 9.658e-002 2 -3.843266909954277e-001 2.112e-001 2.735e-002 3 -2.529596248184762e-001 1.314e-001 7.965e-003 4 -1.677278487209687e-001 8.523e-002 2.343e-003 5 -1.115548118641202e-001 5.617e-002 6.920e-004 6 -7.429254153014103e-002 3.726e-002 2.047e-004 7 -4.950555324952729e-002 2.479e-002 6.063e-005 8 -3.299695812958521e-002 1.651e-002 1.796e-005 9 -2.199597566746803e-002 1.100e-002 5.320e-00610 -1.466339248533952e-002 7.333e-003 1.576e-00611 -9.775419823703882e-003 4.888e-003 4.671e-00712 -6.516894652077655e-003 3.259e-003 1.384e-00713 -4.344581058363196e-003 2.172e-003 4.100e-00814 -2.896382816341328e-003 1.448e-003 1.215e-00815 -1.930920527663641e-003 9.655e-004 3.600e-009>>

xvect ci mostra la lenta convergenza.

Punto b: metodo di bisezione

a=-1;b=1;nmax=15;toll=1e-12;fun='sin(x)-0.5*sin(2*x)';[xvect,xdiff,fx,it,p,c]=bisezione(a,b,nmax,toll,fun); Superato il numero massimo di iterazioni Numero di Iterazioni : 15 Radice calcolata : -6.1035156250000000e-005 Ordine stimato : 1.0000000000000000 Fattore di riduzione : 0.5000000000000000

Tabella riassuntiva - bisezioneiter=1:it;fprintf('%2d %23.15e %15.3e %15.3e\n',[iter' xvect xdiff fx]') 1 0.000000000000000e+000 5.000e-001 0.000e+000 2 -5.000000000000000e-001 2.500e-001 5.869e-002 3 -2.500000000000000e-001 1.250e-001 7.691e-003 4 -1.250000000000000e-001 6.250e-002 9.728e-004 5 -6.250000000000000e-002 3.125e-002 1.220e-004 6 -3.125000000000000e-002 1.563e-002 1.526e-005 7 -1.562500000000000e-002 7.813e-003 1.907e-006 8 -7.812500000000000e-003 3.906e-003 2.384e-007 9 -3.906250000000000e-003 1.953e-003 2.980e-00810 -1.953125000000000e-003 9.766e-004 3.725e-00911 -9.765625000000000e-004 4.883e-004 4.657e-01012 -4.882812500000000e-004 2.441e-004 5.821e-01113 -2.441406250000000e-004 1.221e-004 7.276e-01214 -1.220703125000000e-004 6.104e-005 9.095e-01315 -6.103515625000000e-005 3.052e-005 1.137e-013>>

Confronto tra i due metodi

Il valore ottenuto con la bisezione è più vicino alla soluzione vera x = 0!!!

Newton

Bisezione

Iter xvect xdiff fx....15 -1.930920527663641e-003 9.655e-004 3.600e-009

Iter xvect xdiff fx....15 -6.103515625000000e-005 3.052e-005 1.137e-013

Motivo di tale comportamento?

BisezioneNumero di Iterazioni : 15 Radice calcolata : -6.1035156250000000e-005 Ordine stimato : 1.0000000000000000 Fattore di riduzione : 0.5000000000000000

NewtonNumero di Iterazioni : 15 Radice calcolata : -1.9309205276636407e-003 Ordine stimato : 0.9999897626825431 Fattore di riduzione : 0.6666198373211456

La radice è multipla?

sin 0.5sin 2 0 0

cos cos 2 0 0

sin 2sin 2 0 0

cos 4cos 2 0 3 0

f x x x f

f x x x f

f x x x f

f x x x f

x = 0 è radice multipla con molteplicità m = 3.

Grafico di f ( x )

fplot('sin(x)-0.5*sin(2*x)’,[-1,1])gridtitle('Andamento di f(x)= sin(x)- 0.5*sin(2*x)')

Anche il graficomostra che la radice x=0 è multipla

Modifica metodo di Newton: m noto

x0=-1; nmax=15; toll=1e-12;fun='sin(x)-0.5*sin(2*x)';dfun='cos(x)-cos(2*x)';mol=3;[xvect,xdiff,fx,it,p,c]=newton_m(x0,nmax,toll,fun,dfun,mol);Arresto per azzeramento di dfun

)(

)(1

k

kkk xf

xfmxx

Iter xvect xdiff fx....3 7.266362938868759e-010 1.634e-003 0.000e+000

Modifica del problema:

( )0,

( )

f xF x

f x

fun='(sin(x)-0.5*sin(2*x))./(cos(x)-cos(2*x))';

dfun='1-((sin(x)-0.5*sin(2*x)).*(-sin(x)+2*… sin(2*x)))./(cos(x)-cos(2*x)).^2';

0f x

F xf x

, ,f x x a b

1 1 2

( ) ( )* = -

( ) ( ) ( )

k k kk k k k

k k k k

F x f x f xx x x x

F x f x f x f x

L’applicazione del metodo di Newton al problema implica:

Metodo di Newton (*): risultati

[xvect,xdiff,fx,it,p,c]=newton(x0,nmax,toll,fun,dfun);Warning: Divide by zero.> In C:\analisi_numerica\prog_matlab_new\eq_non_lin\NEWTON.M at line 70Warning: Divide by zero.> In C:\analisi_numerica\prog_matlab_new\eq_non_lin\NEWTON.M at line 62 Numero di Iterazioni : 5 Radice calcolata : NaN Iter xvect xdiff fx 0 -1.000000000000000e+000 0.000e+000 4.044e-001 1 -3.108031246559614e-001 6.892e-001 1.053e-001 2 -9.943189700371891e-003 3.009e-001 3.314e-003 3 -3.276824228314046e-007 9.943e-003 1.092e-007 4 -1.029850733179902e-010 3.276e-007 NaN 5 NaN NaN NaN

Semplificazione della F ( x )

2

sin 1 cos ; 1 cos 1 2cos

sin 1 cos( ) per 0

( ) 1 cos 1 2cos

cos 2 , 0 0

1 2cos

f x x x f x x x

x xf xF x x

f x x x

xF x F

x

Poiché: sin 2 2sin cosx x x

fplot('sin(x)/(1+2*cos(x))’,[-1,1])gridtitle(‘F(x)=sin(x)/(1+2*cos(x))’)

Grafico di F(x)

fplot('(cos(x)+2)/(1+2*cos(x))^2’,[-1,1])gridtitle(‘DF(x)=(cos(x)+2)/(1+2*cos(x))^2’)

Grafico di F’(x)

fplot('(7+2*cos(x)) *sin(x)/(1+2*cos(x))^3',[-1,1])gridtitle('F''''(x)= (7+2*cos(x)) *sin(x)/(1+2*cos(x))^3)')

La derivata seconda non ha segno costante per la convergenza del metodo di Newton occorre prendere vicino alla soluzione!!

Grafico di F’’(x)

0x

Risultati del problema modificatometodo Newton (*)

x0=0.7;fun='sin(x)./(1+2*cos(x))';dfun='(cos(x)+2)./(1+2*cos(x)).^2';[xvect,xdiff,fx,it,p,c]=newton(x0,nmax,toll,fun,dfun);

Numero di Iterazioni : 5 Radice calcolata : 0.0000000000000000e+000 Ordine stimato : 3.0020776453233839 Fattore di riduzione : 0.3387128554185560

x0=-0.8;[xvect,xdiff,fx,it,p,c]=newton(x0,nmax,toll,fun,dfun);

Numero di Iterazioni : 5 Radice calcolata : 0.0000000000000000e+000 Ordine stimato : 3.0052937675224145 Fattore di riduzione : 0.3450700538122051

Utilizzo della function ‘fzero’ di Matlab

>> alfa=fzero('(x^2-4)*cos(x)+4*x*sin(x)',0.8)

alfa =

8.052563483762232e-001

>> options=optimset('Tolx',1e-15);

>> x=fzero('sin(x)-0.5*sin(2*x)',-0.1,options)

x =

1.050493220693351e-008

21. 4 cos 4 sin 0f x x x x x

2. s i n 0.5 s i n 2 0f x x x

Il risultato è meno preciso!!

Esercizio 2: Sistema non lineareDato il sistema non lineare :

2 2

2

4 00 :

0

x y yF x

y x

a) separare graficamente le soluzioni;b) utilizzare il metodo di punto fisso per approssimare la soluzione situata nel primo quadrante;c) applicare quindi il metodo di Newton per risolvere lo stesso problema.d) Confrontare i risultati ottenuti nei due casi.

Punto a: separazione grafica0422 yyx

02 xyCirconferenza:C=(0,2),R = 2

Parabola simm.rispetto l’asse x

La soluzione D = [1,2] x [1,2]

fplot('2-sqrt(4-x^2)', [0,2]),gridhold onfplot('2+sqrt(4-x^2)', [0,2])fplot('sqrt(x)',[0,4]), hold off%y=0:.01:4; %x2=0:.01:4;%x1=sqrt(-%y.^2+4*y);%y2=sqrt(x2);%plot(x1,y,x2,y2),grid

Grafico di F ( x ) in [1,2] x [1,2]fplot('2-sqrt(4-x^2)', [1,2]),grid

hold on

fplot('2+sqrt(4-x^2)', [1,2])

fplot('sqrt(x)',[1,2])axis([1 2 1 2])hold off

Punto b: metodo di punto fisso

)(

)(4

),(

),(

2

12

2

1

xgxy

ygyyx

yxg

yxgxGx

1. continua e derivabile

2. ( )

G

G D D

Convergenza del metodo di punto fisso

1 1

1 2

2 2

,3. 0 1 :

,

g g

x yg gJ x D

x y g g

x y

,Indichiamo: maxi j i j

x y DJ

Verifica delle ipotesi: Hp.1, Hp.2

2 1, 2 [1,2]g D

Hp.1 Le funzioni

1

2

2Hp.2 0

4

dg yy D

dy y y

1 3,2 1,2g D

g1 crescente

g2 è anche crescente

11 2,g g C D

Verifica delle ipotesi: Hp. 3

2

12 212

20

4 2, 0, 0

1 40

2

y

y y yJ J J

y y

x

12 12 12 123/ 22

40 1 1 / 3 1

4J J J

y y

12 21

1max , 1

3

3 / 221 21 21 21

10 1 1 / 2 1

4J J J

x

Il metodo di punto fissoconverge!!

Istruzioni metodo punto fissox0=[1 1];nmax=30; toll=1.e-6;fun=strvcat('x(1)^2+x(2)^2-4*x(2)','-x(1)+x(2)^2'); % x,y sono in x(1),x(2)g=strvcat('sqrt(4*x(2)-x(2)^2)','sqrt(x(1))');[xvect,xdiff,fx,it]=Punto_fissoxs(x0,nmax,toll,fun,g);

....while (it<nmax) & (norm(res_x,inf)>=toll) xap=x; for k=1:n x_new(k)=eval(g(k,:));% x(k)=x_new(k); % da aggiungere per avere la % soluzione con metodo in serie end x=xap;....

Risultati e tabella: metodo in parallelo

Numero di iterazioni : 15 Radice calcolata: 1.9010802796691053e+000 1.3787965146159946e+000 iter=0:it;tab=[iter' xvect xdiff fx];fprintf('%2d %19.15f %19.15f %13.3e %13.3e\n',tab')it soluzione xdiff fx

0 1.000000000000000 1.000000000000000 0.000e+000 2.000e+000 1 1.732050807568877 1.000000000000000 7.321e-001 7.321e-001 2 1.732050807568877 1.316074012952492 3.161e-001 5.322e-001.... .... .... .... ....14 1.901079828717214 1.378796514615995 1.380e-006 1.715e-00615 1.901080279669105 1.378796514615995 4.510e-007 4.510e-007

Risultati e tabella: metodo in serie

it soluzione xdiff fx

0 1.000000000000000 1.000000000000000 0.000e+000 2.000e+000 1 1.732050807568877 1.316074012952492 7.321e-001 5.322e-001 2 1.879426839289333 1.370921893941932 1.474e-001 7.202e-002 3 1.898489066726062 1.377856693102030 1.906e-002 8.677e-003 4 1.900772923229405 1.378685215424248 2.284e-003 1.030e-003 5 1.901043907559104 1.378783488282009 2.710e-004 1.221e-004 6 1.901076023090326 1.378795134561450 3.212e-005 1.447e-005 7 1.901079828717214 1.378796514615995 3.806e-006 1.715e-006 8 1.901080279669105 1.378796678146965 4.510e-007 2.032e-007

Numero di iterazioni : 8 % sono di meno!

Radice calcolata:

1.9010802796691053e+000

1.3787966781469649e+000

Punto c: metodo di Newton

Indichiamo:

2 1 2

0

,( ), det 0 ,

,

Th. : ( )

f fF C D J x y D

x y

I D x I

1k k kJ x h F x

kkk yxx ,

: 0,F D

1 1 perk k kx x h k

Hp.

La soluzione del sistema:

Istruzioni metodo di Newtonx0=[1 1];nmax=30; toll=1.e-6;fun=strvcat('x(1)^2+x(2)^2-4*x(2)','-x(1)+x(2)^2');Jac=strvcat('2*x(1)','2*x(2)-4','-1','2*x(2)');

% Jacobiana in forma di vettore di stringhe

[xvect,xdiff,fx,it]=newtonxs(x0,nmax,toll,fun,Jac);

Numero di iterazioni : 6 % Sono in numero < 8!!Radice calcolata: 1.9010803402881398e+000 1.3787967001295518e+000

Tabella riassuntiva del metodo di Newton

it soluzione xdiff fx

0 1.000000000000000 1.000000000000000 0.000e+000 2.000e+000 1 3.000000000000000 2.000000000000000 2.000e+000 5.000e+000 2 2.166666666666667 1.541666666666667 8.333e-001 9.045e-001 3 1.927083333333333 1.395833333333333 2.396e-001 7.867e-002 4 1.901399523614772 1.379015003483898 2.568e-002 9.425e-004 5 1.901080391293387 1.378796735902042 3.191e-004 1.495e-007 6 1.901080340288140 1.378796700129552 5.101e-008 2.665e-015

iter=0:it;

tab=[iter' xvect xdiff fx];

fprintf('%2d %19.15f %19.15f %13.3e %13.3e\n’, tab')