a
6
lezione - laboratorio
Corso di Laurea ING. MECCANICA
a.a 2004-2005
La scatola dei colori
Avevo una scatola di colori
Ogni colore splendeva di gioia
Avevo una scatola di colori
Caldi e delicati e vivi.
Non avevo il rosso per le ferite e per il sangue,
Non avevo il nero per l’orfano,
Non avevo il bianco per il viso dei morti,
Non avevo il giallo per le sabbie infuocate,
Avevo l’arancione per la gioia e per la vita,
Avevo il verde per germogli e fioriture,
Avevo il blu per limpidi cieli azzurri,
Avevo il rosa per i sogni e per il riposo.
Mi sedetti e dipinsi
La PACE
Tali Shurek
Esercizio 1
Data l’equazione:
f  x   sin  x   0.5sin  2 x   0, x   1, 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.
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
Radice calcolata
Ordine stimato
Fattore di riduzione
:
:
:
:
15
-1.9309205276636407e-003
0.9999897626825431
0.6666198373211456
Tabella riassuntiva - Newton
iter=0:it;
fprintf('%2d
0
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
>>
%23.15e
%15.3e %15.3e\n',[iter' xvect xdiff fx]')
-1.000000000000000e+000
-5.955642027988763e-001
-3.843266909954277e-001
-2.529596248184762e-001
-1.677278487209687e-001
-1.115548118641202e-001
-7.429254153014103e-002
-4.950555324952729e-002
-3.299695812958521e-002
-2.199597566746803e-002
-1.466339248533952e-002
-9.775419823703882e-003
-6.516894652077655e-003
-4.344581058363196e-003
-2.896382816341328e-003
-1.930920527663641e-003
0.000e+000
4.044e-001
2.112e-001
1.314e-001
8.523e-002
5.617e-002
3.726e-002
2.479e-002
1.651e-002
1.100e-002
7.333e-003
4.888e-003
3.259e-003
2.172e-003
1.448e-003
9.655e-004
3.868e-001
9.658e-002
2.735e-002
7.965e-003
2.343e-003
6.920e-004
2.047e-004
6.063e-005
1.796e-005
5.320e-006
1.576e-006
4.671e-007
1.384e-007
4.100e-008
1.215e-008
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
Radice calcolata
Ordine stimato
Fattore di riduzione
:
15
: -6.1035156250000000e-005
: 1.0000000000000000
: 0.5000000000000000
Tabella riassuntiva - bisezione
iter=1:it;
fprintf('%2d
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
>>
%23.15e
%15.3e %15.3e\n',[iter' xvect xdiff fx]')
0.000000000000000e+000
-5.000000000000000e-001
-2.500000000000000e-001
-1.250000000000000e-001
-6.250000000000000e-002
-3.125000000000000e-002
-1.562500000000000e-002
-7.812500000000000e-003
-3.906250000000000e-003
-1.953125000000000e-003
-9.765625000000000e-004
-4.882812500000000e-004
-2.441406250000000e-004
-1.220703125000000e-004
-6.103515625000000e-005
5.000e-001
2.500e-001
1.250e-001
6.250e-002
3.125e-002
1.563e-002
7.813e-003
3.906e-003
1.953e-003
9.766e-004
4.883e-004
2.441e-004
1.221e-004
6.104e-005
3.052e-005
0.000e+000
5.869e-002
7.691e-003
9.728e-004
1.220e-004
1.526e-005
1.907e-006
2.384e-007
2.980e-008
3.725e-009
4.657e-010
5.821e-011
7.276e-012
9.095e-013
1.137e-013
Confronto tra i due metodi
Newton
Iter
xvect
....
15 -1.930920527663641e-003
xdiff
fx
9.655e-004
3.600e-009
xdiff
fx
3.052e-005
1.137e-013
Bisezione
Iter
xvect
....
15 -6.103515625000000e-005
Il valore ottenuto con la bisezione è più
vicino alla soluzione vera x = 0!!!
Motivo di tale comportamento?
Newton
Numero di Iterazioni
Radice calcolata
Ordine stimato
Fattore di riduzione
:
:
:
:
15
-1.9309205276636407e-003
0.9999897626825431
0.6666198373211456
Bisezione
Numero di Iterazioni
Radice calcolata
Ordine stimato
Fattore di riduzione
:
15
: -6.1035156250000000e-005
: 1.0000000000000000
: 0.5000000000000000
La radice è multipla?
f  x   sin  x   0.5sin  2 x 
 f  0  0
f   x   cos  x   cos  2 x 
 f   0  0
f   x   sin  x   2sin  2 x 
 f   0   0
f   x    cos  x   4cos  2 x   f   0   3  0
x = 0 è radice multipla con molteplicità
m = 3.
Grafico di f ( x )
Anche il grafico
mostra che la radice
x=0 è multipla
fplot('sin(x)-0.5*sin(2*x)’,[-1,1])
grid
title('Andamento di f(x)= sin(x)- 0.5*sin(2*x)')
Modifica metodo di Newton: m noto

f ( xk ) 
 x k 1  x k  m

f ( x k ) 

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
Iter
xvect
....
3
7.266362938868759e-010
xdiff
1.634e-003
fx
0.000e+000
Modifica del problema:
 f   x  ,  x  a, b
f  x
F  x 
0
f  x
f ( x)
 0,
L’applicazione del metodo di Newton al problema F  x  
f ( x )
implica:
 *


F  xk  

 xk 1 =xk 
F   xk  





f ( xk ) f ( xk )


 xk 1  xk 

2
 f ( xk )  f ( xk ) f ( xk ) 


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';
…
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 70
Warning: 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
0 -1.000000000000000e+000
1 -3.108031246559614e-001
2 -9.943189700371891e-003
3 -3.276824228314046e-007
4 -1.029850733179902e-010
5
NaN
xdiff
0.000e+000
6.892e-001
3.009e-001
9.943e-003
3.276e-007
NaN
fx
4.044e-001
1.053e-001
3.314e-003
1.092e-007
NaN
NaN
Semplificazione della F ( x )
Poiché:
sin2 x  2sin x cos x
f  x   sin x 1  cos x  ; f   x   1  cos x 1  2cos x 
f ( x)
F  x 

f ( x )
F  x 
sin x 1  cos x 
1  cos x 1  2cos x 
cos x  2
1  2cos x 
2
, F   0  0
per x  0
Grafico di F(x)
fplot('sin(x)/(1+2*cos(x))’,[-1,1])
grid
title(‘F(x)=sin(x)/(1+2*cos(x))’)
Grafico di F’(x)
fplot('(cos(x)+2)/(1+2*cos(x))^2’,[-1,1])
grid
title(‘DF(x)=(cos(x)+2)/(1+2*cos(x))^2’)
Grafico di F’’(x)
La derivata seconda non ha
segno costante per la
convergenza del metodo di
Newton occorre prendere x 0
vicino alla soluzione!!
fplot('(7+2*cos(x)) *sin(x)/(1+2*cos(x))^3',[-1,1])
grid
title('F''''(x)= (7+2*cos(x)) *sin(x)/(1+2*cos(x))^3)')
Risultati del problema modificato
metodo 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
Radice calcolata
Ordine stimato
Fattore di riduzione
:
:
:
:
5
0.0000000000000000e+000
3.0020776453233839
0.3387128554185560
x0=-0.8;
[xvect,xdiff,fx,it,p,c]=newton(x0,nmax,toll,fun,dfun);
Numero di Iterazioni
Radice calcolata
Ordine stimato
Fattore di riduzione
:
:
:
:
5
0.0000000000000000e+000
3.0052937675224145
0.3450700538122051
1.
Utilizzo della function ‘fzero’ di
Matlab
2
f  x    x  4  cos  x   4 x sin  x   0
>> alfa=fzero('(x^2-4)*cos(x)+4*x*sin(x)',0.8)
alfa =
8.052563483762232e-001
2. f  x   sin  x   0.5 sin  2 x   0
>> options=optimset('Tolx',1e-15);
>> x=fzero('sin(x)-0.5*sin(2*x)',-0.1,options)
x =
1.050493220693351e-008
Il risultato è meno preciso!!
Esercizio 2: Sistema non lineare
Dato il sistema non lineare :
 x 2  y 2  4 y  0
F  x  0 :  2
 y  x  0
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 grafica
x 2  y 2  4 y  0 Circonferenza:C=(0,2),R =
2
y  x  0 Parabola simm.rispetto l’asse x
fplot('2-sqrt(4-x^2)',
[0,2]),grid
hold on
fplot('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
La soluzione
  D = [1,2] x [1,2]
2
Grafico di F ( x ) in [1,2] x [1,2]
fplot('2-sqrt(4x^2)',
[1,2]),grid
hold on
fplot('2+sqrt(4x^2)', [1,2])
fplot('sqrt(x)',[
1,2])
axis([1 2 1 2])
hold off
Punto b: metodo di punto fisso
2
 g1 ( x, y ) 
x

4
y

y
 g1 ( y )

x  Gx   

 g 2 ( x, y) 
 y  x  g 2 ( x)
Convergenza del metodo di punto fisso
1. G continua e derivabile
 g1
2. G ( D )  D
 x
 g1 , g2 
3.  0    1 : J 
  
 g2
 x, y 
 x

g1 
y 
 
g2 

y 
Indichiamo:  i j  max J i j
 x , y  D
x  D
Verifica delle ipotesi: Hp.1, Hp.2
Hp.1 Le funzioni g1 , g2  C
1
 D
dg1
2 y
Hp.2

 0 y  D  g1 crescente
2
dy
4y  y
g1  D    3, 2  1, 2


g2 è anche crescente  g2
 D   1,
2   [1, 2]

Verifica delle ipotesi: Hp. 3
2 y

 0
J 
 1

2 x
 
J12

2 
4y  y 
2 y
, J12 
 0, J 21  0

4 y  y2
0


4
4 y  y 
 
J 21
2
1
4x
3/ 2
   max
3/ 2
 0  J12  J12 1  12  1/ 3  1
 0  J 21  J 21 1   21  1/ 2  1

12 ,  21  
1
3
 1 Il metodo di punto fisso
converge!!
Istruzioni metodo punto fisso
x0=[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
0
1
2
....
14
15
soluzione
1.000000000000000
1.732050807568877
1.732050807568877
....
1.901079828717214
1.901080279669105
1.000000000000000
1.000000000000000
1.316074012952492
....
1.378796514615995
1.378796514615995
xdiff
0.000e+000
7.321e-001
3.161e-001
....
1.380e-006
4.510e-007
fx
2.000e+000
7.321e-001
5.322e-001
....
1.715e-006
4.510e-007
Risultati e tabella: metodo in serie
Numero di iterazioni :
8 % sono di meno!
Radice calcolata:
1.9010802796691053e+000
1.3787966781469649e+000
it
0
1
2
3
4
5
6
7
8
soluzione
1.000000000000000
1.732050807568877
1.879426839289333
1.898489066726062
1.900772923229405
1.901043907559104
1.901076023090326
1.901079828717214
1.901080279669105
1.000000000000000
1.316074012952492
1.370921893941932
1.377856693102030
1.378685215424248
1.378783488282009
1.378795134561450
1.378796514615995
1.378796678146965
xdiff
0.000e+000
7.321e-001
1.474e-001
1.906e-002
2.284e-003
2.710e-004
3.212e-005
3.806e-006
4.510e-007
fx
2.000e+000
5.322e-001
7.202e-002
8.677e-003
1.030e-003
1.221e-004
1.447e-005
1.715e-006
2.032e-007
Punto c: metodo di Newton
Indichiamo:
x
k 



k
k
 x ,y

 : F    0,   D
Hp.
  f1 , f 2  
F  C (D ), det  J 
 0 
  x, y  
2
Th.  I    D :  x
 0
La soluzione del sistema: J
x h
x
 k 1
 x
k
 k 1
h
 x, y   D
 I ( )
k
 
 k 1
 F
x 
per k  
k 
Istruzioni metodo di Newton
x0=[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
iter=0:it;
tab=[iter' xvect xdiff fx];
fprintf('%2d
tab')
it
0
1
2
3
4
5
6
%19.15f
%19.15f %13.3e %13.3e\n’,
soluzione
1.000000000000000
3.000000000000000
2.166666666666667
1.927083333333333
1.901399523614772
1.901080391293387
1.901080340288140
1.000000000000000
2.000000000000000
1.541666666666667
1.395833333333333
1.379015003483898
1.378796735902042
1.378796700129552
xdiff
0.000e+000
2.000e+000
8.333e-001
2.396e-001
2.568e-002
3.191e-004
5.101e-008
fx
2.000e+000
5.000e+000
9.045e-001
7.867e-002
9.425e-004
1.495e-007
2.665e-015