1

7a lezione - laboratorio

a.a 2004-2005

Corso di Laurea ING. MECCANICA

2

Si utilizzi il metodo di Runge Kutta 4 per calcolare la soluzione del problema del pendolo:

0

0 0 0 0

sin , tmaxt g t t t

t t

Esercizio 1: Problema del Pendolo

Problema di Cauchy del 2° ordine

nell’incognita t

g: accelerazione di gravità

3

Problema equivalente del 1° ordine

0 0

0 0

sin*

x t y tt x t

y t g x tt x t y t

x tt y t

y t

Risolvendo (*) in si ottengono: x (t) e quindi (t), valore dell’angolo istante per istante, e , velocità. y t t

Sostituzioni Problema equivalente

0 , tmaxt

4

Si assume:

Si analizzano i risultati che si ottengono applicando il metodo di Runge Kutta 4.

Dati del problema e metodo di soluzione

5.00

00

0

0

y

x

0 0 tmax 10

9.81260 valore unificato

100

t

g

n

5

Istruzioni: Runge Kutta 4 t0=0;tmax=10;n=100;y0=[0 0.5];f1='y(2)';f2='-9.81260*sin(y(1))';f=strvcat(f1,f2);[T,Y]=Rungekutta4(t0,tmax,n,y0,f);x=Y(:,1);y=Y(:,2);plot(T,x)axis([0 10 min(x) max(x)])title('Valori istantanei dell''angolo-RK4')xlabel('Tempo (s)')ylabel('Angolo (rad)')tabella

6

7

File “tabella.m”% Il file permette di scrivere una tabella in cui si riportano i valori

% in uscita da un problema differenziale di Cauchy di tipo vettoriale.

% Se si vuole riportare anche l'errore, vanno aggiunte altre colonne

clc

a=[T,x,y];

s='----------------------------------------------';

disp(s)

fprintf('tempo\t\t\t angolo\t\t\t\t velocità\n')

disp(s)

n=length(T);

for i=1:10:n

fprintf('%6.3f %16.6e %16.6e \n',a(i,1:3))

end

8

Risultati del Metodo RK4

---------------------------------------------tempo angolovelocità--------------------------------------------- 0.000 0.000000e+000 5.000000e-001 1.000 2.289632e-003 -4.999165e-001 2.000 -4.578396e-003 4.997301e-001 3.000 6.865820e-003 -4.994409e-001 4.000 -9.151435e-003 4.990491e-001 5.000 1.143477e-002 -4.985546e-001 6.000 -1.371536e-002 4.979577e-001 7.000 1.599273e-002 -4.972584e-001 8.000 -1.826642e-002 4.964570e-001 9.000 2.053595e-002 -4.955537e-001 10.000 -2.280087e-002 4.945486e-001

9

Movimento del punto P

Il punto P ha coordinate:

ttsinP cos,

Per costruire la rappresentazione del punto P in movimento, determiniamo le coordinate di P, cioè: (sin ( x ( i )), - cos ( x ( i ))), i = 0,1,…,n .

x

y

P

N.B. l’asse y è orientato verso il basso, cioè in modo contrario della normale rappresentazione.

10

Metodo di Runge Kutta 4: simulazione

………………………y0=[0 3];[T,Y]=Rungekutta4(t0,tmax,n,y0,f);theta=Y(:,1);n=length(T);plot(0,-1,'or');for i=1:n x(i)=sin(theta(i)); y(i)=-cos(theta(i)); plot(x(i),y(i),'ob',[0,x(i)],[0,y(i)],'r'); axis([-1 1 -1.5 .5]) pause(.25)endtitle('Oscillazioni del pendolo – RK4')xlabel('Ascissa del punto P')ylabel('Ordinata del punto P')

11

12

Si determini la traiettoria ed il tempo di attraversamento di un fiume largo 2 km e con velocità della corrente di modulo s, da parte di un battello che si muove con velocità relativa (rispetto all’acqua) di modulo v, e che, partendo da un punto a valle (o a monte) del punto di attracco, si dirige sempre verso tale punto. Si utilizzi il metodo di Eulero.

Si esaminerà il caso di partenza da un punto a valle. Lo studente può studiare l’altro caso.

Esercizio 2: Moto del Battello

13

ssyy

dx/dtdx/dt

VV

vv

dy/dtdy/dt

xxAA

PPssyy

dx/dtdx/dt

VV

vv

dy/dtdy/dt

xxAA

PP

2 2 kmkm

0.1

0.1

kmkm

14

1.00

20

22

22

ysyx

vy

dt

dy

xyx

vx

dt

dx

Modello del problema

dx dyV i j

dt dt Velocità assoluta battello:

15

Eulero scalare:

Metodo di Eulero

Eulero vettoriale:

Eulero vettoriale applicato in serie:

iiii ztfhzz ,*1

[T,Y]=Eulero(t0,tmax,n,y0,f)

f=strvcat(f1,f2);[T,Y]=Eulero(t0,tmax,n,y0,f)

iiiii

iiiii

yxtfhyy

yxtfhxx

,,*

,,*

21

11

5.0;0

*

*

,,*

,,*

00

11

1

121

11

yx

xsinghyy

yhxx

yxtfhyy

yxtfhxxiii

iii

iiiii

iiiii

[T,X,Y]=Eulero2(t0,tmax,n,x0,y0,f1,f2);

16

Soluzione del problema e simulazione del moto

t0=0;tmax=0.5; %primo valore di tentativoy0=[2 0.1];

f=strvcat('-5*y(1)/sqrt(y(1)^2+y(2)^2)‘, ...'-5*y(2)/sqrt(y(1)^2+y(2)^2)+3');n=50;[T,Y]=Eulero(t0,tmax,n,y0,f);x=Y(:,1); y=Y(:,2);figure(1)axis([-0.1 2 -0.2 0.6])hold onfor i=1:length(T) plot(x(i),y(i),'*b') %*=simbolo,b=blu pause(0.30)endhold offtitle([' Traiettoria del battello - tmax = ' num2str(tmax)])xlabel(' X [km]');ylabel(' Y [km]')grid

17

L’estremo finale viene calcolato per tentativi.

18

19

Risultati con tmax=0.65 T X Y 0.000 2.000000e+000 1.000000e-001 0.013 1.935081e+000 1.357541e-001 0.026 1.870240e+000 1.702052e-001 0.039 1.805508e+000 2.033141e-001 0.052 1.740916e+000 2.350406e-001 ... ... ... 0.559 1.919500e-002 1.821415e-001 0.572 1.238269e-002 1.564994e-001 0.585 7.255724e-003 1.307019e-001 0.598 3.652893e-003 1.048019e-001 0.611 1.388678e-003 7.884131e-002 0.624 2.439726e-004 5.285139e-002 0.637 -5.607722e-005 2.685208e-002 0.650 7.966688e-005 8.522218e-004

20

Sistema differenziale non lineare di ordine 1

0 0

1

2

0 0

P t k P t c P t Q t

Q t k Q t d P t Q t

P P Q Q

Esercizio 3: Problema Preda-Predatore

P(t): prede, Q(t): predatori. Si trovi la soluzione del problema per assumendo: k1 = 2; k2 = 10; c = 0.001; d = 0.002; P0=5000; Q0=100.

0,3t

21

Soluzione analitica del problema

2 22

1 1 1

k d P Q k d PdQ k Q d PQ Q

dP k P c PQ k c Q P P k c Q

1 2k c Q k d PdQ dP

Q P

1ln 2ln cost.

1ln 2ln cost.

k Q c Q k P d P

k Q c Q k P d P

Integrando i due membri si ottiene l’Int. Gen.

cost. si calcola imponendo le condizioni iniziali.

Variabili separate

22

Istruzioni: Metodo Eulero esplicito

e grafici di P( t ), Q( t )t0=0;tmax=3;y0=[5000 100]; % P,Q sono in y(1), y(2)f=strvcat('2*y(1)- 0.001*y(1)*y(2)',...'-10*y(2)+ 0.002*y(1)*y(2)');str1='Eulero';n=300; h=tmax./n;[T,Y]=Eulero(t0,tmax,n,y0,f);plot(T,Y(:,1),’b’,T,Y(:,2),’g’); gridtitle(['Risultati del metodo di ',str1])legend('prede','predatori')title(['Risultati del metodo di ',str1,...' con h =' num2str(h)])xlabel('Tempo')

23

Grafici: Metodo di Eulero Esplicito

24

Istruzioni per grafico nel piano P Qt0=0;tmax=3;

y0=[5000 100]; % P,Q sono in y(1), y(2)

f=strvcat('2*y(1)- 0.001*y(1)*y(2)','-10*y(2)+ .002*y(1)*y(2)');

str1='Eulero';N=[60 300]; h=tmax./N;

n=N(1);

[T,Y]=Eulero(t0,tmax,n,y0,f);

plot(Y(:,1),Y(:,2),'r')

hold on

n=N(2);

[T,Y]=Eulero(t0,tmax,n,y0,f);

plot(Y(:,1),Y(:,2),'b',10/0.002,2/0.001,'*'),grid

title(['Risultati del metodo di ',str1])

xlabel('Prede'),ylabel('Predatori');grid

legend(['h1 = ' num2str(h(1))],['h2 = 'num2str(h(2))])

25

Rappresentazione nel piano PQ

Il punto segnatocon * è il puntocritico che si ottiene ponendo:

1 0

2 0

P t k P cPQ

Q t k Q dPQ

*

*

2 /

1/

P k d

Q k c

26

Istruzioni del Metodo Heun Grafici di P( t ), Q( t )

t0=0;tmax=3;y0=[5000 100];f=strvcat('2*y(1)- 0.001*y(1)*y(2)',...'-10*y(2)+ 0.002*y(1)*y(2)');str1=‘Heun';n=300; h=tmax./n; [T,Y]=Heun(t0,tmax,n,y0,f);plot(T,Y(:,1),T,Y(:,2)); gridlegend('prede','predatori')title(['Risultati del metodo di ',str1,... ' con h =' num2str(h)])xlabel('Tempo')

27

Grafici prede, predatoriMetodo di Heun

28

Soluzione nel piano PQ:simulazione movimento

t0=0;tmax=3;y0=[5000 100];

f=strvcat('2*y(1)- 0.001*y(1)*y(2)','-10*y(2)+ 0.002*y(1)*y(2)');

n=300; h=tmax./n;

[T,Y]=Heun(t0,tmax,n,y0,f);

plot(10/0.002,2/0.001,’*r’,y0(1),y0(2),’*b’)

hold on

for i=1:n

plot(Y(i,1),Y(i,2),’*b’)

pause(0.25)

end

title('Simulazione nel piano PQ - Heun')

29

Simulazione: metodo di Heun

30

1 1

1 1

10 1 0 1 0

1

Y t Y tx t x t y t

y t x t y t

x y Y

Esercizio 4: Sistema differenziale lineare del 1° ordine

a) Stabilire se il problema ammette soluzione unica e se è ben condizionato.

b)Calcolare la soluzione con il metodo di Eulero Implicito, n1=100 e n2=200;

n1, n2 = numero sottointervalli.

0,10t

31

sint cost

sint cost

t

t

x t e

y t e

c) Si calcoli l’errore nei nodi sapendo che la soluzione vera è:

d) Si confronti la soluzione vera calcolata nei nodi al punto c), con quella approssimata ottenuta applicando il metodo di Runge-Kutta 4 con n1 = 100.

it

Esercizio 4: quesiti c,d

32

a: esistenza, unicità della soluzione

1

2

, ,,

, ,

x t x t y t f t x t y tY t F t Y

y t x t y t f t x t y t

1 1,

1 1FJ Y

2FJ Y

Esiste unica la soluzione del problema.

33

a: stabilità e condizionamento

La matrice dei coefficienti è i suoi autovalori si trovano risolvendo l’equazione caratteristica:

1 1

1 1A

1,2 1,21 1 1 0 1 0i e

il sistema è asintoticamente stabile e quindi ben condizionato.

34

n=input('n = '); % numero sottointervalli t0=0;tmax=10;coeff=[-1 1 0; -1 -1 0];x0=-1;y0=1;[T,X,Y]=Eulsis(t0,tmax,n,x0,y0,coeff);str1='Eulero Implicito';

b: metodo di Eulero Implicito

disp(['Risultati del metodo di ',str1]);t=T;Xv=exp(-t).*(sin(t)-cos(t));Yv=exp(-t).*(sin(t)+cos(t));tabella2

Stampa dei risultati

35

c: file “tabella2.m”% Il file scrive una tabella in cui si riportano i valori in uscita da

% un problema differenziale di Cauchy di tipo vettoriale.

% Si vuole riportare l'errore nella X e nella Y, quindi vanno

% considerate le colonne relative.

%

ErrX= abs(X-Xv);ErrY= abs(Y-Yv);a=[T,X,Y, ErrX,ErrY];s='------------------------------------------';disp(s)fprintf(' T X Y ErrX ErrY \n');disp(s)n=length(T);for i=1:10:n %stampa ogni 10 valori fprintf('%6.3f %16.6e %16.6e %10.2e %10.2e \n',a(i,:))end

36

----------------------------------------------------------------------------------------------------- T X Y ErrX ErrY----------------------------------------------------------------------------------------------------- 0.000 -1.000000e+000 1.000000e+000 0.00e+000 0.00e+000 1.000 6.326408e-002 5.194194e-001 4.75e-002 1.11e-002 2.000 1.657577e-001 1.000365e-001 1.36e-002 3.33e-002 3.000 6.695044e-002 -2.547604e-002 1.06e-002 1.68e-002 4.000 7.847666e-003 -2.531598e-002 9.74e-003 5.17e-004 5.000 -5.585724e-003 -8.060362e-003 2.79e-003 3.51e-003 6.000 -3.622299e-003 -2.110340e-004 5.50e-004 1.90e-003 7.000 -8.876485e-004 1.007195e-003 7.99e-004 2.79e-004 8.000 9.098508e-005 4.883277e-004 2.90e-004 2.05e-004 9.000 1.630219e-004 8.486890e-005 2.80e-007 1.46e-004 10.000 6.190752e-005 -2.813839e-005 4.85e-005 3.47e-005

Tabella Eulero Implicito: n1 = 100

37

% Grafico soluzione vera ed approssimataplot(T,Xv,T,X);grid xlabel('Tempo (s)')title(’Metodo di Eulero implicito - n=100')legend('Xvera','Xapp')

38

Errore per la x(t): Eulero Implicito

39

Tabella Eulero Implicito: n2 = 200-------------------------------------------------------------- T X Y ErrX ErrY-------------------------------------------------------------- 0.000 -1.000000e+000 1.000000e+000 0.00e+000 0.00e+000 0.500 -2.615168e-001 8.176226e-001 2.00e-002 5.44e-003 1.000 8.623577e-002 5.138799e-001 2.46e-002 5.55e-003 1.500 1.894160e-001 2.532959e-001 1.74e-002 1.49e-002 2.000 1.726328e-001 8.400313e-002 6.75e-003 1.73e-002 2.500 1.165047e-001 -2.675497e-003 1.62e-003 1.40e-002 3.000 6.211849e-002 -3.383809e-002 5.80e-003 8.42e-003 3.500 2.410848e-002 -3.553024e-002 6.42e-003 3.34e-003 4.000 3.128919e-003 -2.587447e-002 5.02e-003 4.12e-005 ...... 7.000 -5.159904e-004 1.197169e-003 4.28e-004 8.94e-005 7.500 5.446362e-005 7.894289e-004 2.73e-004 7.89e-005 8.000 2.488899e-004 4.108081e-004 1.32e-004 1.28e-004 8.500 2.485198e-004 1.524551e-004 3.64e-005 1.12e-004 9.000 1.764843e-004 1.315848e-005 1.32e-005 7.47e-005 9.500 9.888436e-005 -4.197206e-005 2.99e-005 3.83e-005 10.000 4.168455e-005 -5.014193e-005 2.83e-005 1.27e-005

40

Costruzione tabella riassuntivan1=100; % numero sottointervalli

t0=0;tmax=10;coeff=[-1 1 0; -1 -1 0]; x0=-1;y0=1;

[T1,X1,Y1]=Eulsis(t0,tmax,n1,x0,y0,coeff);

t=T1;

Xv=exp(-t).*(sin(t)-cos(t)); Yv=exp(-t).*(sin(t)+cos(t));

ErrX1= abs(X1-Xv);ErrY1= abs(Y1-Yv);

n2=200; % numero sottointervalli [T2,X2,Y2]=Eulsis(t0,tmax,n2,x0,y0,coeff);

t=T2;

Xv=exp(-t).*(sin(t)-cos(t)); Yv=exp(-t).*(sin(t)+cos(t));

ErrX2= abs(X2-Xv);ErrY2= abs(Y2-Yv);

a=[T1,ErrX1,ErrX2(1:2:end),ErrY1,ErrY2(1:2:end)];

fprintf(' T ErrX1 ErrX2 ErrY1 ErrY2\n');

fprintf('%6.3f %10.2e %10.2e %10.2e %10.2e \n',a(1:10:end,:)')

41

Tabella riassuntiva dei risultati T ErrX1 ErrX2 ErrY1 ErrY2

0.000 0.00e+000 0.00e+000 0.00e+000 0.00e+000

1.000 4.75e-002 2.46e-002 1.11e-002 5.55e-003

2.000 1.36e-002 6.75e-003 3.33e-002 1.73e-002

3.000 1.06e-002 5.80e-003 1.68e-002 8.42e-003

4.000 9.74e-003 5.02e-003 5.17e-004 4.12e-005

5.000 2.79e-003 1.28e-003 3.51e-003 1.92e-003

6.000 5.50e-004 3.86e-004 1.90e-003 9.42e-004

7.000 7.99e-004 4.28e-004 2.79e-004 8.94e-005

8.000 2.90e-004 1.32e-004 2.05e-004 1.28e-004

9.000 2.80e-007 1.32e-005 1.46e-004 7.47e-005

10.000 4.85e-005 2.83e-005 3.47e-005 1.27e-005

42

MODIFICHE ai files EULSIS e TRAPEZI

.. …. .. . %file EULSIS.M

while t < tmax

tn=x1+h*tn1;

x1=mat\tn;

X=[X;x1(1)];

Y=[Y;x1(2)];

t=t+h;

if abs(tmax-t)<1.e-13

t=tmax;

end

T=[T;t];

end

.. …. .. . %file TRAPEZI.M

while t<tmax

c1_old=eval(coeff(1,:));

c2_old=eval(coeff(2,:));

t=t+h;

if abs(tmax-t)<1.e-13

t=tmax;

end

c1_new=eval(coeff(1,:));

c2_new=eval(coeff(2,:));

.. …. .. .

end

43

Soluzione con movimento nel piano XY

n=200;t0=0;tmax=10;coeff=[-1 1 0;-1 -1 0];x0=-1;y0=1; [T,X,Y]=Eulsis(t0,tmax,n,x0,y0,coeff);plot(0,0,'or',x0,y0,'*g') % (0,0) punto % di stazionarietàhold onfor i=1:nplot(X(i),Y(i),'ob')pause(0.25)end

44

Grafico nel piano XY

45

n=input('n = '); % numero sottointervalli t0=0;tmax=10;y0=[-1 1]; % x,y sono in y(1), y(2)f=strvcat('-y(1)+y(2)','-y(1)-y(2)'); [T,Y]=RungeKutta4(t0,tmax,n,y0,f);X=Y(:,1); Y=Y(:,2);str1='Runge-Kutta';

d: metodo di Runge-Kutta 4

disp(['Risultati del metodo di ',str1]);t=T;Xv=exp(-t).*(sin(t)-cos(t));Yv=exp(-t).*(sin(t)+cos(t));tabella

Stampa dei risultati

46

Tabella Metodo RK4: n1 = 100

-------------------------------------------------------------- T X Y ErrX ErrY -------------------------------------------------------------- 0.000 -1.000000e+000 1.000000e+000 0.00e+000 0.00e+000 1.000 1.107954e-001 5.083239e-001 1.62e-006 2.12e-006 2.000 1.793787e-001 6.673883e-002 6.65e-007 1.85e-006 3.000 5.631372e-002 -4.226311e-002 1.05e-006 2.42e-007 4.000 -1.889800e-003 -2.583285e-002 3.79e-007 3.71e-007 5.000 -8.372432e-003 -4.549641e-003 4.95e-008 2.39e-007 6.000 -3.072525e-003 1.687461e-003 1.01e-007 3.86e-008 7.000 -8.833836e-005 1.286537e-003 3.73e-008 2.74e-008 8.000 3.807013e-004 2.830635e-004 1.21e-009 1.94e-008 9.000 1.632948e-004 -6.158686e-005 7.03e-009 3.92e-009 10.000 1.339237e-005 -6.279076e-005 2.90e-009 1.55e-009

47

% Grafico soluzione vera ed approssimataplot(T,Xv,T,X);grid xlabel('Tempo (s)')title(’Metodo di Runge-Kutta - n=100')legend('Xvera','Xapp')

48

Errore per la x(t): Metodo RK4