Post on 19-Jan-2016
description
transcript
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')