1 8 a bis lezione di laboratorio Laurea Specialistica in Ingegneria Matematica Ingegneria dei...

Post on 01-May-2015

217 views 0 download

Transcript of 1 8 a bis lezione di laboratorio Laurea Specialistica in Ingegneria Matematica Ingegneria dei...

1

8a bis lezione di laboratorio

Laurea Specialistica in

Ingegneria MatematicaIngegneria dei Sistemi Energetici

Laurea Specialistica in

Ingegneria MatematicaIngegneria dei Sistemi Energetici

a.a. 2007-2008

2

Si utilizzi il metodo RK4 per calcolare la soluzione del problema del pendolo:

2

02

0 0 0 0

sin , tmaxd

g t t td t

dt t

d t

Esercizio 1: Problema di Cauchy

Problema di Cauchy del 2° ordine

nell’incognita t

g: accelerazione di gravità

3

Problema equivalente del 1° ordine

12

1

21 1

2

2 1 0 02

2 0 0

sin*

d yy t

dtt y td yd yd g y ty t d td t dty td ydy td t d t

Risolvendo (*) in si ottengono:

e quindi (t), valore dell’angolo istante per istante, e , velocità. 2

dy

dt

Sostituzioni Problema equivalente

0 , tmaxt 1y t

4

Inizialmente si assume:

e si analizzano i risultati che si ottengono applicando il metodo di Runge-Kutta4.

Dati del problema

5.00

00

20

10

y

y

0 0 tmax 10

9.81260 valore unificato

100

t

g

n

5

Istruzioni Metodo RK4t0=0;tmax=10;n=100;y0=[0 0.5];g=num2str(9.81260, '%9.5f');f1='y(2)';f2=[g,'*(-sin(y(1)))'];f=strvcat(f1,f2);[T,Y]=Rungekutta4(t0,tmax,n,y0,f);x=Y(:,1);y=Y(:,2);plot(T,x)title('Valori istantanei dell''angolo -Rungekutta4')xlabel('Tempo (s)')ylabel('Angolo (rad)')tabella

6

7

File “tabella.m”% Il file scrive 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];

ar=a(1:10:end,:);

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

disp(s)

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

disp(s)

fprintf('%6.3f %16.6e %16.6e \n',ar’)

8

Risultati del Metodo RK4

---------------------------------------------tempo angolo velocità--------------------------------------------- 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

Considerazioni sul problema

• Il sistema è conservativo:

21 2 2 1

1, 1 cos C

2E y y y g y

Energia cinetica Energia potenziale

10

1. ; è il livello minimo di energia e dalla (1) essendo i due termini non negativi, deve essere , e la (1) si riduce al punto di equilibrio (critico) (0,0).

2. ; sulla base di (2), otteniamo, considerando i due segni, una

traiettoria chiusa intorno all’origine. Questi livelli di energia, corrispondono al moto oscillatorio intorno al punto .

3. ; per la (2) abbiamo una curva chiusa nel piano contenente 4 orbite:

i due punti di equilibrio ed altre due orbite che connettono

tali punti, dette eteroclite.

4. ; in questo caso la (2) individua due orbite che non si chiudono. Tali

livelli di energia corrispondono al moto di rivoluzione completa del pendolo che

descrive così l’intera circonferenza.

22 1

11 1 cos 0,

2E y g y C C

2 12 2 cos ,y C g g y

0C 1 20, 0y y

0 2C g

2C g

2C g

,0

11

Piano delle fasi y1y2

-4 -2 0 2 4-3

-2

-1

0

1

2

3

E=2k

E=C3

E=C2E=C1

y1

y2k=g

12

Metodo RK4: piano delle fasi

% Il piano delle fasi e’ stato ottenuto variando% tmax e le condizioni iniziali%t0=0;tmax=100;n=1000;y0=[pi-0.01 0]; g=num2str(9.81260, '%9.5f');f1='y(2)';f2=[g,'*(-sin(y(1)))'];f=strvcat(f1,f2);[T,Y]=Rungekutta4(t0,tmax,n,y0,f);x=Y(:,1);y=Y(:,2);plot(x,y)title('Piano delle fasi-RK4')xlabel('x')ylabel('y')

13

14

Solver ode45 e piano delle fasi

t0=0;tmax=100;y0=[pi-0.01 0];[T,Y]=ode45('pendolo',[t0 tmax],y0);x=Y(:,1);y=Y(:,2);plot(x,y)title('Piano delle fasi-ode45')xlabel('x')ylabel('y')

t0=0;tmax=100;y0=[pi-0.01 0];[T,Y]=ode45('pendolo',[t0 tmax],y0);x=Y(:,1);y=Y(:,2);plot(x,y)title('Piano delle fasi-ode45')xlabel('x')ylabel('y')

function f=pendolo(t,y)f=[y(2)-9.81260*sin(y(1))];

function f=pendolo(t,y)f=[y(2)-9.81260*sin(y(1))];

15

16

Solver ode45 con utilizzo di options e piano delle fasi

t0=0;tmax=100;y0=[pi-0.01 0];options=odeset('RelTol',1e-6,'AbsTol',1e-7);[T,Y]=ode45('pendolo',[t0 tmax],y0,options);x=Y(:,1);y=Y(:,2);plot(x,y)title('Piano delle fasi-ode45')xlabel('x')ylabel('y')

17

18

Movimento del punto PNel sistema fissato in figura, il punto P ha coordinate:

sin ,cost t

Per costruire però la rappresentazione del punto P in movimento utilizzando Matlab, occorre determinare le coordinate di P in tale sistema, quindi: (sin ( x ( i )),-cos ( x ( i ))), i = 0,1,…,n .

x

y

P

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

19

Metodo RK4 - simulazione

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 – Runge-Kutta4') xlabel('Ascissa del punto P')ylabel('Ordinata del punto P')

20

21

Un paracadutista del peso di 80 kg viene lanciato da un aeroplano ad un’altezza di 600 m. Dopo 5 s il paracadute si apre. Nel seguente modello la funzione y rappresenta l’altezza del paracadutista ad ogni istante t:

Esercizio 2

2

2

( )

(0) 600 , (0) 0 / ;

d y tg

d t m

d yy m m s

d t

2

2

( )

(0) 600 , (0) 0 / ;

d y tg

d t m

d yy m m s

d t

g=9.81260 è l’accelerazione di gravità,

22

m=80 è la massa del paracadutista, è la resistenza dell’aria proporzionale alquadrato della velocità con 2 diversi coefficienti di proporzionalità a seconda che si consideri prima o dopo l’apertura del paracadute:

( )t

21

22

( '( )) 5( )

( '( )) 5 .

K y t t st

K y t t s

1 -Si consideri dapprima il caso Sapendo che la soluzione analitica è

1 0K 2 0K

2( ) 600,2

gy t t

23

si chiede: a) a quale altezza si trova il paracadutista dopo 5s?b) quanto tempo impiega il paracadutista per raggiungere il terreno? c) qual’ è la velocità di impatto?

2 -Si consideri poi il caso e si determini per tentativi il tempo necessario per l’impatto a terra utilizzando un metodo multistep.

1 2

1 4,

150 150K K

Dopo aver calcolato la soluzione approssimata,si risponda alle domande a), b) e c) del puntoprecedente.

24

Punto 1: lancio senza attritoclear allclc% Lancio senza attrito % altezza dopo 5sg=9.8126;altezza=600;ta=5;ya=-g*ta^2/2+altezza %tempo di impatto: 0=-g/2*tf^2+altezzatf=sqrt(2*altezza/g)% y'(t) calcolata al tempo finaley1=-g*tf t=0:0.2:tf;y=-g*t.^2/2+altezza;

Altezza dopo 5s = 4.773425000000000e+002tempo di impatto = 11.05855991281188velocità di impatto = -1.085132250004579e+002

Risultati

25

Istruzioni per la simulazione

axis([0 12 0 600])hold onfor i=1:length(t) plot(t(i),y(i),'*') pause(.25)end

hold offtitle('Altezza in funzione del tempo-Assenza di attrito')xlabel('tempo')ylabel('altezza')

26

27

Punto 2: metodo multistepAdams-Bashforth (A-B)

% Lancio in presenza di attritoclear allclcg=9.81260;gs=num2str(g,'%9.5f');m=80;ms=num2str(m);k=[1 4]/150; % da 0 a 5 secondi:t0=0;t5=5;y0=[600 0];ks1=num2str(k(1),'%24.16e');f1='y(2)';f2=['-',gs,'+(',ks1,'*y(2).^2)/',ms];f=strvcat(f1,f2);n=50;h=(t5-t0)/n;[t_inn,y_inn]= Rungekutta4(t0,t0+3*h,3,y0,f);[T1,Y1]=Adamsbf(t0,t5,n,y_inn,f);

% Lancio in presenza di attritoclear allclcg=9.81260;gs=num2str(g,'%9.5f');m=80;ms=num2str(m);k=[1 4]/150; % da 0 a 5 secondi:t0=0;t5=5;y0=[600 0];ks1=num2str(k(1),'%24.16e');f1='y(2)';f2=['-',gs,'+(',ks1,'*y(2).^2)/',ms];f=strvcat(f1,f2);n=50;h=(t5-t0)/n;[t_inn,y_inn]= Rungekutta4(t0,t0+3*h,3,y0,f);[T1,Y1]=Adamsbf(t0,t5,n,y_inn,f);

28

Applicazione di A-B dopo 5s

% dopo 5 secondi: tmax=11.3637643; % è stato calcolato per tentativiy5=Y1(end,:);ks2=num2str(k(2),'%24.16e');f1='y(2)'; f2=['-',gs,'+(',ks2,'*y(2).^2)/',ms];f=strvcat(f1,f2);n=150;h=(tmax-t5)/n;[t_inn,y_inn]= Rungekutta4(t5,t5+3*h,3,y5,f);[T2,Y2]=Adamsbf(t5,tmax,n,y_inn,f);

% dopo 5 secondi: tmax=11.3637643; % è stato calcolato per tentativiy5=Y1(end,:);ks2=num2str(k(2),'%24.16e');f1='y(2)'; f2=['-',gs,'+(',ks2,'*y(2).^2)/',ms];f=strvcat(f1,f2);n=150;h=(tmax-t5)/n;[t_inn,y_inn]= Rungekutta4(t5,t5+3*h,3,y5,f);[T2,Y2]=Adamsbf(t5,tmax,n,y_inn,f);

29

Costruzione delle tabelletab1=[T1 Y1];tab1r=tab1(1:5:end,:);s='_______________________________________________ ';disp(s)fprintf('\n TEMPO ALTEZZA VELOCITA''\n')disp(s)fprintf(' %10.5f %20.10f %20.10f\n',tab1r')fprintf('\n\n')tab2=[T2 Y2];tab2r=tab2(1:10:end,:);disp(s)fprintf('\n TEMPO ALTEZZA VELOCITA''\n')disp(s)fprintf(' %10.5f %20.10f %20.10f\n',tab2r')disp(s)

30

Tabella (da 0s a 5s)_____________________________________________________

TEMPO ALTEZZA VELOCITA’_____________________________________________________ 0.00000 600.0000000000 0.0000000000 0.50000 598.7734667893 -4.9059656963 1.00000 595.0943685169 -9.8099262298 1.50000 588.9642084400 -14.7098797166 2.00000 580.3854892578 -19.6038308142 2.50000 569.3617090292 -24.4897939651 3.00000 555.8973554760 -29.3657966060 3.50000 539.9978986930 -34.2298823317 4.00000 521.6697822891 -39.0801140039 4.50000 500.9204129902 -43.9145767939 5.00000 477.7581487393 -48.7313811492

31

Tabella (da 5s al tempo d’impatto)__________________________________________________ TEMPO ALTEZZA VELOCITA’__________________________________________________ 5.00000 477.7581487393 -48.7313811492 5.42425 456.2757618953 -52.5316215452 5.84850 433.1928931664 -56.2758262841 6.27275 408.5339767395 -59.9609176276 6.69700 382.3247032605 -63.5840510944 7.12125 354.5919197763 -67.1426203313 7.54551 325.3635279901 -70.6342602128 … … … 10.09101 120.2549361827 -90.0765456317 10.51526 81.4056373438 -93.0536164828 10.93951 41.3097005084 -95.9533672841 11.36376 0.0000010018 -98.7754907131__________________________________________________

32

Istruzioni per la simulazioneaxis([0 12 0 600])hold onfor i=1:length(T1) plot(T1(i),Y1(i,1),'*r') pause(.25)endfor i=1:length(T2) plot(T2(i),Y2(i,1),'*b') pause(.25)endhold offtitle('Altezza in funzione del tempo-Adams-Bashforth')xlabel('tempo')ylabel('altezza')

33

34

ode113 Adams-Bashforth-Moultonn=50;t0=0;tmax=5;h=(tmax-t0)/n;global g m kg= 9.81260;m=80;k=1/150; [T1,Y1]=ode113('paracadute',[t0:h:tmax],[600,0],[ ],k);tmax1=11.3637643;n=150;h=(tmax1-T1(end))/n;global kk=4/150;[T2,Y2]=ode113('paracadute',[T1(end):h:tmax1],[Y1(end,:)],[],k);

function f=paracadute(t,y)global g m kf=[y(2)-g+(k*y(2)^2)/m];

function f=paracadute(t,y)global g m kf=[y(2)-g+(k*y(2)^2)/m];

35

Risultati con ode113: da t0=0 a tmax=5s _____________________________________________ TEMPO ALTEZZA VELOCITA’ _____________________________________________

0.00000 600.0000000000 0.0000000000 0.50000 598.7734667988 -4.9059656971 1.00000 595.0943685927 -9.8099262325 1.50000 588.9642086456 -14.7098797214 2.00000 580.3854896221 -19.6038308212 2.50000 569.3617095516 -24.4897939741 3.00000 555.8973561512 -29.3657966169 3.50000 539.9978993290 -34.2298823433 4.00000 521.6697829093 -39.0801140170 4.50000 500.9204135962 -43.9145768085 5.00000 477.7581493297 -48.7313811653

36

Risultati dopo 5s – ode113 _________________________________________________ TEMPO ALTEZZA VELOCITA’ __________________________________________________

5.00000 477.7581493297 -48.7313811653 5.42425 456.2757588114 -52.5316211051 5.84850 433.1928899696 -56.2758258422 6.27275 408.5339646237 -59.9609168796 6.69700 382.3247058183 -63.5840508266 7.12125 354.5919316388 -67.1426203541 7.54551 325.3635384887 -70.6342601703 7.96976 294.6683922545 -74.0568481943 8.39401 262.5361927432 -77.4085041130 8.81826 228.9973906228 -80.6875883812 9.24251 194.0830813906 -83.8926987113 9.66676 157.8249037315 -87.0226653705 10.09101 120.2549451618 -90.0765455475 10.51526 81.4056463446 -93.0536164012 10.93951 41.3097095387 -95.9533672050 11.36376 0.0000100676 -98.7754906355

37

Esercizio 3Si consideri un saltatore in lungo con il punto di partenza del salto situatoall’origine degli assi coordinati. Il problema della traiettoria del salto è:

Le costanti del problema sono

La resistenza aerodinamica è data da ed il coefficiente di resistenza

e l’area della sezione trasversale del saltatore sono rispettivamente: è la densità dell’aria diversa

ad una certa altezza rispetto al livello del mare.1 - Dalla natura del problema si desume che nessuna delle variabili del problema può

crescere indefinitamente. Stabilire quindi se il problema ammette soluzione unica .

0

cos sin0 0 0 0

0 8 0cos sin

d x d yv t t v t t

d t d t x y

D t v vd g dt g t

d t v t d t m

0

cos sin0 0 0 0

0 8 0cos sin

d x d yv t t v t t

d t d t x y

D t v vd g dt g t

d t v t d t m

9.81260, 80 .g m kg 2 2

2

c s d x d yD

d t d t

20.72, 0.50 ,c s m

38

2 - Si considerino i due casi: a)

b)

Si deve determinare per tentativi il valore tmax per cui si ottiene y( tmax )=0.Si costruisca un file MATLAB che, una volta avviato:• permetta di dare in input il numero di sottointervalli, n=100, in [t0, tmax];• calcoli la soluzione approssimata utilizzando il metodo di Runge-Kutta4 ed

un solver di MATLAB (ode23);• faccia visualizzare una tabella riassuntiva che riporti l’intestazione: tmax distanza e riporti su due righe il valore di tmax per cui y(tmax)=0, trovato nei due

casi e la corrispondente distanza ed utilizzi i seguenti formati di stampa: 3 cifre decimali e formato virgola fissa per i valori di tmax; 12 cifre decimali e formato virgola fissa per le distanze di atterraggio .

3 - Mediante subplot con 2 finestre grafiche su una riga, riporti le traiettorie nei 2 casi corredate dell’ indicazione della distanza e della velocità di impatto.

Si commentino i risultati.

300.94 / , 10 / ;kg m v m s

301.29 / , 10 / .kg m v m s

39

Soluzione con RK4: istruzioni clear all;clcc=0.72;s=0.50;g=9.81260;m=80;v0=10;gs=num2str(g,'%10.5f');d1=-c*s/(2*m);t0=0;y0=[0 0 pi/8 v0];n=100;f1='y(4)*cos(y(3))';f2='y(4)*sin(y(3))'; f3=[gs, '/(-y(4))*cos(y(3))'];% caso arho=0.94;tmax1=input('tmax= ');% tmax1=0.7780 determinato per tentativic1=num2str(d1*rho,'%25.16e');f4=[c1, '*y(4)^2-', gs, '*sin(y(3))'];f=strvcat(f1,f2,f3,f4);[T1,Y1]=Rungekutta4(t0,tmax1,n,y0,f);a1=[T1,Y1]; disp('RISULTATI CASO a)')fprintf(' t x y theta ni \n')fprintf('%6.3f %14.6e %14.6e %14.6e %14.6e\n',a1(1:10:end,:)') % caso brho=1.29;tmax2=input('tmax= '); % tmax2=0.7770 determinato per tentativic2=num2str(d1*rho,'%25.16f');f4=[c2, '*y(4)^2-', gs, '*sin(y(3))'];f=strvcat(f1,f2,f3,f4);[T2,Y2]=Rungekutta4(t0,tmax2,n,y0,f);a2=[T2,Y2];

40

disp(' ');disp('RISULTATI CASO b)')fprintf('t x y theta ni \n')fprintf('%6.3f %14.6e %14.6e %14.6e %14.6e\n', a2(1:10:end,:)') % confronto caso a) e caso b)tab=[a1(end,1:2); a2(end,1:2)];fprintf('tmax distanza \n')fprintf('%6.3f %16.12f \n',tab') % Graficisubplot(121), plot(Y1(:,1),Y1(:,2),'*'),grid,axis([0 8 0 0.8])title(['a) distanza=' num2str(Y1(end,1)) ' vel='

num2str(Y1(end,4))]); xlabel('x');ylabel('y')subplot(122), plot(Y2(:,1),Y2(:,2),'o'),grid, axis([0 8 0 0.8])title(['b) distanza=' num2str(Y2(end,1)) 'vel='

num2str(Y2(end,4))]);xlabel('x');ylabel('y')

Soluzione con RK4: istruzioni

41

Risultati caso a)

__________________________________________________________________________ t x y theta ni ________________________________________________________________

0.000 0.000000e+000 0.000000e+000 3.926991e-001 1.000000e+001 0.078 7.181931e-001 2.678043e-001 3.201124e-001 9.717474e+000 0.156 1.435238e+000 4.758335e-001 2.437382e-001 9.489762e+000 0.233 2.151164e+000 6.241899e-001 1.642286e-001 9.320445e+000 0.311 2.865992e+000 7.129700e-001 8.244516e-002 9.212326e+000 0.389 3.579736e+000 7.422664e-001 -5.823280e-004 9.167171e+000 0.467 4.292404e+000 7.121693e-001 -8.372662e-002 9.185523e+000 0.545 5.003996e+000 6.227685e-001 -1.658518e-001 9.266625e+000 0.622 5.714502e+000 4.741554e-001 -2.459041e-001 9.408474e+000 0.700 6.423909e+000 2.664244e-001 -3.229865e-001 9.607998e+000 0.778 7.132194e+000 -3.249795e-004 -3.964061e-001 9.861308e+000

__________________________________________________________________________ t x y theta ni ________________________________________________________________

0.000 0.000000e+000 0.000000e+000 3.926991e-001 1.000000e+001 0.078 7.181931e-001 2.678043e-001 3.201124e-001 9.717474e+000 0.156 1.435238e+000 4.758335e-001 2.437382e-001 9.489762e+000 0.233 2.151164e+000 6.241899e-001 1.642286e-001 9.320445e+000 0.311 2.865992e+000 7.129700e-001 8.244516e-002 9.212326e+000 0.389 3.579736e+000 7.422664e-001 -5.823280e-004 9.167171e+000 0.467 4.292404e+000 7.121693e-001 -8.372662e-002 9.185523e+000 0.545 5.003996e+000 6.227685e-001 -1.658518e-001 9.266625e+000 0.622 5.714502e+000 4.741554e-001 -2.459041e-001 9.408474e+000 0.700 6.423909e+000 2.664244e-001 -3.229865e-001 9.607998e+000 0.778 7.132194e+000 -3.249795e-004 -3.964061e-001 9.861308e+000

42

Risultati caso b)

__________________________________________________________________________ t x y theta ni ________________________________________________________________

0.000 0.000000e+000 0.000000e+000 3.926991e-001 1.000000e+001 0.078 7.170537e-001 2.674145e-001 3.201857e-001 9.711879e+000 0.155 1.432538e+000 4.750665e-001 2.438455e-001 9.478789e+000 0.233 2.146493e+000 6.230961e-001 1.643246e-001 9.304261e+000 0.311 2.858951e+000 7.116356e-001 8.248136e-002 9.191057e+000 0.389 3.569933e+000 7.408113e-001 -6.546264e-004 9.140909e+000 0.466 4.279450e+000 7.107464e-001 -8.395242e-002 9.154320e+000 0.544 4.987503e+000 6.215633e-001 -1.662691e-001 9.230495e+000 0.622 5.694083e+000 4.733861e-001 -2.465415e-001 9.367386e+000 0.699 6.399170e+000 2.663433e-001 -3.238622e-001 9.561867e+000 0.777 7.102738e+000 5.694180e-004 -3.975284e-001 9.809983e+000

__________________________________________________________________________ t x y theta ni ________________________________________________________________

0.000 0.000000e+000 0.000000e+000 3.926991e-001 1.000000e+001 0.078 7.170537e-001 2.674145e-001 3.201857e-001 9.711879e+000 0.155 1.432538e+000 4.750665e-001 2.438455e-001 9.478789e+000 0.233 2.146493e+000 6.230961e-001 1.643246e-001 9.304261e+000 0.311 2.858951e+000 7.116356e-001 8.248136e-002 9.191057e+000 0.389 3.569933e+000 7.408113e-001 -6.546264e-004 9.140909e+000 0.466 4.279450e+000 7.107464e-001 -8.395242e-002 9.154320e+000 0.544 4.987503e+000 6.215633e-001 -1.662691e-001 9.230495e+000 0.622 5.694083e+000 4.733861e-001 -2.465415e-001 9.367386e+000 0.699 6.399170e+000 2.663433e-001 -3.238622e-001 9.561867e+000 0.777 7.102738e+000 5.694180e-004 -3.975284e-001 9.809983e+000

43

Risultati di Runge-kutta4

44

Soluzione caso a) con il solver ode23 function yp=saltatore(t,y)global g m c s rhod= c*rho*s/2;yp=[y(4)*cos(y(3));y(4)*sin(y(3));-g/y(4)*cos(y(3)); -d/m*y(4)^2-g*sin(y(3))];

clear all; clcglobal m c g s rhom=80;c=0.72;s=0.5;rho=0.94;g = 9.81260;tmax=0.778;[T,Y]=ode23('saltatore',[0,tmax],[0 0 pi/8 10]);

Al tempo finale,t= 0.7780 ,ottenuto per tentativi, si ottengono i seguenti valori delle componenti: 7.1322 -0.0005 -0.3964 9.8615

45

Risultati ode23 – caso a)

__________________________________________________________________________ t x y theta ni ________________________________________________________________

0 0 0 3.9270e-001 1.0000e+001 8.6591e-006 8.0000e-005 3.3137e-005 3.9269e-001 1.0000e+001 5.1955e-005 4.8000e-004 1.9881e-004 3.9265e-001 9.9998e+000 2.6843e-004 2.4800e-003 1.0269e-003 3.9246e-001 9.9989e+000 1.3508e-003 1.2480e-002 5.1604e-003 3.9147e-001 9.9946e+000 6.7628e-003 6.2476e-002 2.5654e-002 3.8655e-001 9.9734e+000 3.3823e-002 3.1237e-001 1.2378e-001 3.6164e-001 9.8707e+000 1.1162e-001 1.0301e+000 3.6557e-001 2.8734e-001 9.6115e+000 1.8942e-001 1.7466e+000 5.4763e-001 2.0951e-001 9.4088e+000 2.6722e-001 2.4620e+000 6.7005e-001 1.2889e-001 9.2658e+000 3.4502e-001 3.1764e+000 7.3294e-001 4.6434e-002 9.1850e+000 4.2282e-001 3.8897e+000 7.3639e-001 -3.6785e-002 9.1674e+000 5.0033e-001 4.5992e+000 6.8080e-001 -1.1931e-001 9.2130e+000 5.7813e-001 5.3103e+000 5.6585e-001 -2.0066e-001 9.3206e+000 6.5593e-001 6.0204e+000 3.9173e-001 -2.7953e-001 9.4877e+000 7.3373e-001 6.7293e+000 1.5853e-001 -3.5510e-001 9.7110e+000 7.7800e-001 7.1322e+000 -4.8995e-004 -3.9640e-001 9.8615e+000

__________________________________________________________________________ t x y theta ni ________________________________________________________________

0 0 0 3.9270e-001 1.0000e+001 8.6591e-006 8.0000e-005 3.3137e-005 3.9269e-001 1.0000e+001 5.1955e-005 4.8000e-004 1.9881e-004 3.9265e-001 9.9998e+000 2.6843e-004 2.4800e-003 1.0269e-003 3.9246e-001 9.9989e+000 1.3508e-003 1.2480e-002 5.1604e-003 3.9147e-001 9.9946e+000 6.7628e-003 6.2476e-002 2.5654e-002 3.8655e-001 9.9734e+000 3.3823e-002 3.1237e-001 1.2378e-001 3.6164e-001 9.8707e+000 1.1162e-001 1.0301e+000 3.6557e-001 2.8734e-001 9.6115e+000 1.8942e-001 1.7466e+000 5.4763e-001 2.0951e-001 9.4088e+000 2.6722e-001 2.4620e+000 6.7005e-001 1.2889e-001 9.2658e+000 3.4502e-001 3.1764e+000 7.3294e-001 4.6434e-002 9.1850e+000 4.2282e-001 3.8897e+000 7.3639e-001 -3.6785e-002 9.1674e+000 5.0033e-001 4.5992e+000 6.8080e-001 -1.1931e-001 9.2130e+000 5.7813e-001 5.3103e+000 5.6585e-001 -2.0066e-001 9.3206e+000 6.5593e-001 6.0204e+000 3.9173e-001 -2.7953e-001 9.4877e+000 7.3373e-001 6.7293e+000 1.5853e-001 -3.5510e-001 9.7110e+000 7.7800e-001 7.1322e+000 -4.8995e-004 -3.9640e-001 9.8615e+000

46

Esercizio 4: Problema d’esame Appello 22/09/2005

Si consideri un satellite in orbita ellittica attorno ad un pianeta di massa M. Supponiamo che il pianeta sia posto nell’origine delriferimento cartesiano bidimensionale e il problema sianormalizzato in modo che GM=1, con G costante di gravitazioneuniversale.Il moto del satellite è allora regolato dalle equazioni:

Sia T il periodo di rivoluzione del satellite.1 - Si stabilisca se il problema ammette soluzione unica.2 - Per la terza legge di Keplero risulta dove a è il semiasse maggiore dell’orbita.

22 2 3/ 2

2

22 2 3/ 2

2

/( )

/( ) .

d xx x y

d t

d yy x y

d t

22 2 3/ 2

2

22 2 3/ 2

2

/( )

/( ) .

d xx x y

d t

d yy x y

d t

322 4 aT

47

Quesiti 2a, 2bSi costruisca un file MATLAB: Cognome_studente_matricola.m che, una volta avviato:

a) faccia visualizzare una schermata con i dati personali ed una breve presentazione del problema;

b) risolva il problema con il metodo RK4 ed utilizzando ODE45 nei casi

Nel primo caso l’orbita è circolare centro l’origine ed a = 1; nel secondo caso è ancora circolare con a = 2, infine nel terzo caso è un’ellisse con centro (-1,0) e semiasse maggiore a = 2.

(0) 1, (0) 0, (0) 0, (0) 1;d x d x

x yd t d t

(0) 2, (0) 0, (0) 0, (0) 1/ 2;d x d y

x yd t d t

(0) 1, (0) 0, (0) 0, (0) 6 / 2;d x d y

x yd t d t

48

Quesiti 2c, 2d, 3

c) faccia visualizzare tre tabelle riassuntive, una per ogni caso, che riportino ogni 10 nodi:

Intestazione: tempo soluzione derivata_soluzione; con i seguenti formati di stampa:

3 cifre decimali e formato virgola fissa per i valori dei nodi; 8 cifre decimali e virgola fissa per la soluzione e la derivata ;

d) mediante il comando subplot, con 6 finestre grafiche, si ritrovino per via grafica (corredata di label e titolo) le affermazioni del punto b) e si tracci, nei tre casi, l’andamento della soluzione in funzione del tempo.

3 - Si commentino i risultati.

49

Considerazioni teoriche

Il campo gravitazionale è conservativo, pertanto l’energia totale del satellite si conserva:

2

2 2

2 2

1

21

0 0 02 0 0

tot

tot

mE t mv GM C

dm

C E m x y GMx y

v = velocità del pianeta, m = massa del satellite,d = distanza satellite-pianeta

50

Caso in esamePer le condizioni assegnate nei tre casi, tenendo conto che

GM=1, abbiamo:

(1)

(2)

(3)

10 0

2 1 21 1

0 1 02 2 2 2 2

6 30 1 0

2 4 1 4

tot

tot

tot

m mE m

m mE m

mE m m

Pertanto in ognuno dei tre casi si hanno orbite circolari/ellittiche.

51

2a) file XXXX_XXX.m

clear alldisp('Cognome e nome studente: XXXX_XXX')disp('Numero di matricola: XXXX')disp(' Corso di Laurea: XXXX ')disp(' ')disp(' Questo programma consente di calcolare ')disp(' e visualizzare la soluzione ')disp(' del seguente sistema: ')disp(' ')disp(' / x''''(t)=-x/(x^2+y^2)^(3/2)' )disp(' \ y''''(t)=-y/(x^2+y^2)^(3/2)')disp(' con diverse condizioni iniziali')disp(' ')disp(' utilizzando Runge-Kutta4 e il solver ode45. ')

52

Terzo caso-Metodo RK4t0=0;tmax=2*pi*sqrt(2^3);f1='y(2)';f2='-y(1)/(y(1)^2+y(3)^2)^(3/2)';f3='y(4)';f4='-y(3)/(y(1)^2+y(3)^2)^(3/2)';f=strvcat(f1,f2,f3,f4);y0=[1 0 0 sqrt(6)/2]; n=100;[T,Y]=Rungekutta4(t0,tmax,n,y0,f);X=Y(:,1); X1=Y(:,2); YV=Y(:,3);YV1=Y(:,4);str1='Runge-Kutta';disp(['Risultati del metodo di ',str1]);tab=[T X YV X1 YV1];tab_10=tab(1:10:end,:);s='------------------------------------------------------------------'disp(s)fprintf(' tempo soluzione derivata_soluzione \

n')disp(s)fprintf(' %7.3f %12.8f %12.8f %12.8f %12.8f \n',tab_10')subplot(2,1,1),plot(T,X,T,YV),grid,title('soluzioni x(t) e y(t)')xlabel('t'),ylabel('x,y')subplot(2,1,2),plot(X,YV),grid,title('orbita')xlabel('x'),ylabel('y')

53

Risultati RK4 –terzo caso

__________________________________________________________________ tempo soluzione derivata_soluzione __________________________________________________________________

0.000 1.00000000 0.00000000 0.00000000 1.22474487 1.777 -0.03265426 1.51593990 -0.81632406 0.39064849 3.554 -1.35404660 1.70460966 -0.63935333 -0.09962284 5.331 -2.28448065 1.32746219 -0.41023654 -0.29773434 7.109 -2.82363255 0.71084577 -0.19934779 -0.38356134 8.886 -2.99985645 -0.00027980 0.00006125 -0.40826680 10.663 -2.82341431 -0.71137187 0.19947120 -0.38352275 12.440 -2.28404148 -1.32787730 0.41036158 -0.29764433 14.217 -1.35338881 -1.70479154 0.63947005 -0.09943831 15.994 -0.03186270 -1.51561630 0.81630315 0.39107078 17.772 1.00000057 0.00109146 -0.00092295 1.22473755

__________________________________________________________________ tempo soluzione derivata_soluzione __________________________________________________________________

0.000 1.00000000 0.00000000 0.00000000 1.22474487 1.777 -0.03265426 1.51593990 -0.81632406 0.39064849 3.554 -1.35404660 1.70460966 -0.63935333 -0.09962284 5.331 -2.28448065 1.32746219 -0.41023654 -0.29773434 7.109 -2.82363255 0.71084577 -0.19934779 -0.38356134 8.886 -2.99985645 -0.00027980 0.00006125 -0.40826680 10.663 -2.82341431 -0.71137187 0.19947120 -0.38352275 12.440 -2.28404148 -1.32787730 0.41036158 -0.29764433 14.217 -1.35338881 -1.70479154 0.63947005 -0.09943831 15.994 -0.03186270 -1.51561630 0.81630315 0.39107078 17.772 1.00000057 0.00109146 -0.00092295 1.22473755

54

Traiettoria del satellite nel terzo caso – ode45

[t,y]=ode45('satellite1',[0,2*pi*sqrt(2^3)],[1 0 0 sqrt(6)/2]);plot(y(:,1),y(:,3))title('traiettoria del satellite nel terzo caso')

55

Traiettoria del satellite nel terzo caso – ode45

t0=0;tmax=2*pi*sqrt(2^3);n=100;h=(tmax-t0)/n;[T,Y]=ode45('satellite1',[t0:h:tmax],[1 0 0 sqrt(6)/2]);X=Y(:,1); X1=Y(:,2); YV=Y(:,3);YV1=Y(:,4);tab=[T X YV X1 YV1];

function yp=satellite1(t,y)yp=[y(2) -y(1)/(y(1)^2+y(3)^2)^(3/2) y(4) -y(3)/(y(1)^2+y(3)^2)^(3/2)];

function yp=satellite1(t,y)yp=[y(2) -y(1)/(y(1)^2+y(3)^2)^(3/2) y(4) -y(3)/(y(1)^2+y(3)^2)^(3/2)];

56

Risultati ode45 __________________________________________________________________ tempo soluzione derivata_soluzione __________________________________________________________________

0 1.0000e+000 0 0 1.2247e+000 3.6917e-004 1.0000e+000 4.5214e-004 -3.6917e-004 1.2247e+000 3.0354e-003 1.0000e+000 3.7176e-003 -3.0354e-003 1.2247e+000 2.0468e-002 9.9979e-001 2.5067e-002 -2.0465e-002 1.2245e+000 1.2814e-001 9.9182e-001 1.5652e-001 -1.2727e-001 1.2148e+000 8.2667e-001 6.9837e-001 9.1453e-001 -6.4884e-001 9.0379e-001 1.8860e+000 -1.2192e-001 1.5554e+000 -8.1451e-001 3.4382e-001 3.3389e+000 -1.2145e+000 1.7198e+000 -6.6780e-001 -6.3382e-002 5.1762e+000 -2.2206e+000 1.3680e+000 -4.2939e-001 -2.8747e-001 7.3977e+000 -2.8769e+000 5.9156e-001 -1.6581e-001 -3.9199e-001 9.6191e+000 -2.9687e+000 -3.0719e-001 8.2468e-002 -4.0438e-001 1.1841e+001 -2.5052e+000 -1.1458e+000 3.3780e-001 -3.3481e-001 1.4062e+001 -1.4480e+000 -1.6936e+000 6.1854e-001 -1.2306e-001 1.6078e+001 3.6484e-002 -1.4867e+000 8.1386e-001 4.2689e-001 1.7150e+001 8.2033e-001 -7.2472e-001 5.3819e-001 1.0183e+000 1.7705e+001 9.9977e-001 -9.2192e-002 7.3374e-002 1.2187e+000 1.7772e+001 1.0025e+000 -1.1213e-002 7.5563e-003 1.2221e+000

>

__________________________________________________________________ tempo soluzione derivata_soluzione __________________________________________________________________

0 1.0000e+000 0 0 1.2247e+000 3.6917e-004 1.0000e+000 4.5214e-004 -3.6917e-004 1.2247e+000 3.0354e-003 1.0000e+000 3.7176e-003 -3.0354e-003 1.2247e+000 2.0468e-002 9.9979e-001 2.5067e-002 -2.0465e-002 1.2245e+000 1.2814e-001 9.9182e-001 1.5652e-001 -1.2727e-001 1.2148e+000 8.2667e-001 6.9837e-001 9.1453e-001 -6.4884e-001 9.0379e-001 1.8860e+000 -1.2192e-001 1.5554e+000 -8.1451e-001 3.4382e-001 3.3389e+000 -1.2145e+000 1.7198e+000 -6.6780e-001 -6.3382e-002 5.1762e+000 -2.2206e+000 1.3680e+000 -4.2939e-001 -2.8747e-001 7.3977e+000 -2.8769e+000 5.9156e-001 -1.6581e-001 -3.9199e-001 9.6191e+000 -2.9687e+000 -3.0719e-001 8.2468e-002 -4.0438e-001 1.1841e+001 -2.5052e+000 -1.1458e+000 3.3780e-001 -3.3481e-001 1.4062e+001 -1.4480e+000 -1.6936e+000 6.1854e-001 -1.2306e-001 1.6078e+001 3.6484e-002 -1.4867e+000 8.1386e-001 4.2689e-001 1.7150e+001 8.2033e-001 -7.2472e-001 5.3819e-001 1.0183e+000 1.7705e+001 9.9977e-001 -9.2192e-002 7.3374e-002 1.2187e+000 1.7772e+001 1.0025e+000 -1.1213e-002 7.5563e-003 1.2221e+000

>

57

Risultati ode45

__________________________________________________________________ tempo soluzione derivata_soluzione __________________________________________________________________ 0.000 1.00000000 0.00000000 0.00000000 1.22474487 1.777 -0.03305215 1.51545895 -0.81681670 0.38986649 3.554 -1.35524433 1.70202928 -0.63967240 -0.10100981 5.331 -2.28575099 1.32245141 -0.41004068 -0.29902988 7.109 -2.82417987 0.70382777 -0.19877651 -0.38450148 8.886 -2.99918600 -0.00854207 0.00082137 -0.40871305 10.663 -2.82138094 -0.71991871 0.20019742 -0.38338889 12.440 -2.28097785 -1.33563370 0.41071603 -0.29690160 14.217 -1.35049883 -1.71079189 0.63878331 -0.09836109 15.994 -0.03153727 -1.52088170 0.81392305 0.38992698 17.772 1.00245582 -0.01121327 0.00755633 1.22210410

__________________________________________________________________ tempo soluzione derivata_soluzione __________________________________________________________________ 0.000 1.00000000 0.00000000 0.00000000 1.22474487 1.777 -0.03305215 1.51545895 -0.81681670 0.38986649 3.554 -1.35524433 1.70202928 -0.63967240 -0.10100981 5.331 -2.28575099 1.32245141 -0.41004068 -0.29902988 7.109 -2.82417987 0.70382777 -0.19877651 -0.38450148 8.886 -2.99918600 -0.00854207 0.00082137 -0.40871305 10.663 -2.82138094 -0.71991871 0.20019742 -0.38338889 12.440 -2.28097785 -1.33563370 0.41071603 -0.29690160 14.217 -1.35049883 -1.71079189 0.63878331 -0.09836109 15.994 -0.03153727 -1.52088170 0.81392305 0.38992698 17.772 1.00245582 -0.01121327 0.00755633 1.22210410

58

Grafici terzo caso