77aa - 8 - 8aa lezione di laboratorio lezione di laboratorio

Laurea Specialisticain

Ingegneria Matematica Ingegneria dei Sistemi Energetici

Laurea Specialisticain

Ingegneria Matematica Ingegneria dei Sistemi Energetici

a.a. 2007-2008

Esercizio 1Esercizio 1Prova pratica del 04-09-2007Prova pratica del 04-09-2007

.20',20

],2,0()(')('1)('')(')(1 2

yy

xxyxyxyxyxy

.20',20

],2,0()(')('1)('')(')(1 2

yy

xxyxyxyxyxy

Si consideri il seguente problema di Cauchy:

1- Si verifichi che la funzione è soluzione del problema proposto e si dica, motivando opportunamente la risposta, se tale soluzione è unica, sapendo che y(x) e la sua derivata sono funzioni limitate e strettamente positive nell’intervallo [0, 2].

xexy 2

Quesito 2Quesito 2

2- Si costruisca un file MATLAB: Cognome_studente_matricola.m che, una volta avviatoa) faccia visualizzare una schermata con i dati personali ed una breve presentazione del problemab) determini la soluzione approssimata utilizzando il metodo di Heun ed il metodo di Runge Kutta del quarto ordine, con il passo h = 0.04; c) valuti l’errore relativo nei nodi nei due casi; d) faccia visualizzare una tabella riassuntiva in cui si riporti:

intestazione: x solH solRK errH errRKin cui x contiene i nodi xi, presi ogni 5, solH e solRK

sono le soluzioni approssimate in tali nodi, calcolate con i due metodi ed errH e errRK i corrispondenti errori relativi. Si utilizzino i seguenti formati di stampa:3 cifre decimali e formato virgola fissa per i valori dei nodi,10 cifre decimali e formato virgola fissa per le soluzioni nei due casi,2 cifre decimali e formato floating point per gli errori.

Quesiti 3 e 4Quesiti 3 e 4

3- Si esegua una figura con 2 finestre grafiche; nella prima finestra grafica si riporti la soluzione vera (color blu) e quella approssimata ottenuta con il metodo di Heun (color rosso); nella seconda finestra grafica si riporti la soluzione vera (color blu) e quella approssimata ottenuta con il metodo di Runge Kutta (color verde). Si corredino i grafici di titolo e label.

4- Si commentino i risultati e si specifichi se essi soddisfano la aspettative teoriche.

Risoluzione quesito 2:Risoluzione quesito 2:variabili di input delle functionsvariabili di input delle functionsclear allclct0=0;tmax=2;h=0.04;n=round((tmax-t0)/h); y0=[2 2];f1='y(2)';f2='(1+y(2).^2).*y(2)./(1+y(1).*y(2))';f=strvcat(f1,f2);yveras='2*exp(T)';

Risoluzione quesito 2:Risoluzione quesito 2:applicazione dei metodi e stampeapplicazione dei metodi e stampe

[T,Y1]=Heun(t0,tmax,n,y0,f);[T,Y2]=Rungekutta4(t0,tmax,n,y0,f); yvera=eval(yveras);err1=abs(yvera-Y1(:,1))./abs(yvera);err2=abs(yvera-Y2(:,1))./abs(yvera);tab=[T Y1(:,1) Y2(:,1) err1 err2];tabr=tab(1:5:end,:);% stampa della tabella e graficifprintf(' nodi solH solRK errH errRK \n') fprintf('%7.3f %14.10f %14.10f %10.2e %10.2e\n',tabr')subplot(2,1,1),plot(T,yvera,T,Y1(:,1),'r')xlabel('t'),ylabel('yvera -Heun')title('Soluzione vera - Heun h=0.04')subplot(2,1,2),plot(T,yvera,T,Y2(:,1),'g')xlabel('t'),ylabel('yvera - RK4')title('Soluzione vera - RK4 h=0.04')

Risultati quesito 2: tabellaRisultati quesito 2: tabella

nodi solH solRK errH errRK

0.000 2.0000000000 2.0000000000 0.00e+000 0.00e+000 0.200 2.4426790826 2.4428055062 5.18e-005 4.13e-009 0.400 2.9833405503 2.9836493707 1.04e-004 8.25e-009 0.600 3.6436717793 3.6442375557 1.55e-004 1.24e-008 0.800 4.4501604196 4.4510817835 2.07e-004 1.65e-008 1.000 5.4351568857 5.4365635447 2.59e-004 2.06e-008 1.200 6.6381720177 6.6402336811 3.11e-004 2.48e-008 1.400 8.1074619672 8.1103996994 3.62e-004 2.89e-008 1.600 9.9019638802 9.9060645217 4.14e-004 3.30e-008 1.800 12.0936600235 12.0992944794 4.66e-004 3.71e-008 2.000 14.7704651859 14.7781115880 5.17e-004 4.13e-008

nodi solH solRK errH errRK

0.000 2.0000000000 2.0000000000 0.00e+000 0.00e+000 0.200 2.4426790826 2.4428055062 5.18e-005 4.13e-009 0.400 2.9833405503 2.9836493707 1.04e-004 8.25e-009 0.600 3.6436717793 3.6442375557 1.55e-004 1.24e-008 0.800 4.4501604196 4.4510817835 2.07e-004 1.65e-008 1.000 5.4351568857 5.4365635447 2.59e-004 2.06e-008 1.200 6.6381720177 6.6402336811 3.11e-004 2.48e-008 1.400 8.1074619672 8.1103996994 3.62e-004 2.89e-008 1.600 9.9019638802 9.9060645217 4.14e-004 3.30e-008 1.800 12.0936600235 12.0992944794 4.66e-004 3.71e-008 2.000 14.7704651859 14.7781115880 5.17e-004 4.13e-008

Risultati quesito 3: graficiRisultati quesito 3: grafici

Esercizio 2Esercizio 2Prova pratica del 30-11-2006Prova pratica del 30-11-2006

.2

1)0(',0)0(

,0)0('',1)0(',0)0(

'4'''

]1,0(,'''''' 2

zz

yyy

yzzxz

xzyyyxy

.2

1)0(',0)0(

,0)0('',1)0(',0)0(

'4'''

]1,0(,'''''' 2

zz

yyy

yzzxz

xzyyyxy

Si consideri il problema differenziale di Cauchy

1- Sapendo che nell’intervallo [0,1] le componenti della soluzione e le rispettive derivate fino all’ordine 3 si mantengono limitate, si dica, motivando la risposta, se il problema è ben posto.

Quesito 2 Quesito 2

2- Si 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 utilizzando il metodo Runge Kutta del quarto ordine con i passi h1=0.1, h2=0.05 e h3=0.01;

c) faccia visualizzare una tabella riassuntiva che riporti ogni 5 nodi: intestazione: nodi sol1 sol2 sol3

dove nodi sono i nodi comuni delle partizioni, sol1, sol2, sol3 sono le soluzioni numeriche ottenute con il metodo e valutate nei nodi comuni, con i seguenti formati di stampa:

3 cifre decimali e formato virgola fissa per i valori dei nodi; 4 cifre decimali e formato virgola mobile per le soluzioni.

Quesito 3, 4 e 5Quesito 3, 4 e 5

3- In una stessa figura corredata di label e titolo, si riporti la componente y(x) di ognuna delle tre soluzioni numeriche ottenute con il metodo RK4.

4- Si utilizzi l’ode suite per risolvere il problema proposto, scegliendo ode45 e riportando nella stessa figura del punto precedente, il grafico della componente y(x) della soluzione ottenuta con ode45 ( color rosso e punto).

5- Si commentino i risultati.

Risoluzione quesito 2Risoluzione quesito 2clear allclc t0=0; tmax=1;y0=[0 1 0 0 0.5];h=[0.1 0.05 0.01];n=round((tmax-t0)./h);f1='y(2)';f2='y(3)';f3='t.^2.*y(1).*y(3)-y(1).*y(5)';f4='y(5)';f5='t.*y(4).*y(5)+4*y(2)';f=strvcat(f1,f2,f3,f4,f5); tab=[];s='%7.3f';for i=1:length(n) step=round(h(1)/h(i)); [T,Y]=Rungekutta4(t0,tmax,n(i),y0,f); tab=[tab Y(1:step:end,[1,4])]; subplot(1,2,1),plot(T,Y(:,1));hold on subplot(1,2,2),plot(T,Y(:,4));hold on s=[s,' %12.4e %12.4e'];endtab=[T(1:step:end) tab];tabr=tab(1:5:end,:);

Risoluzione quesiti 3 e 4Risoluzione quesiti 3 e 4

fprintf(' t \t\t\t\t\t sol1_h1 \t\t\t\t sol2_h2 \t\t\t\t sol3_h3 \n')fprintf([s,' \n'],tabr')[T,Y] = ode45('fun5',[t0 tmax],y0);tab=[T Y(:,:)];tab=tab(1:4:end,:);fprintf(' t solode45 \n')fprintf('%7.3f %12.4e %12.4e %12.4e %12.4e %12.4e \n', tab')subplot(1,2,1),plot(T,Y(:,1),'r.')xlabel('t');ylabel('y(t)');title('Componente y(t)')subplot(1,2,2),plot(T,Y(:,4),'r.')xlabel('t');ylabel('z(t)');title('Componente z(t)')

function f5=fun5(t,y)f5=[y(2) y(3) t.^2.*y(1).*y(3)-y(1).*y(5) y(5) t.*y(4).*y(5)+4*y(2)];

Risultati quesito 2: tabellaRisultati quesito 2: tabella

t sol1_h1 sol2_h2 sol3_h3

0.000 0.0000e+000 0.0000e+000 0.0000e+000 0.0000e+000 0.0000e+000 0.0000e+000 0.500 4.9659e-001 7.5910e-001 4.9659e-001 7.5910e-001 4.9659e-001 7.5910e-001 1.000 9.0431e-001 2.9880e+000 9.0430e-001 2.9882e+000 9.0430e-001 2.9882e+000

t sol1_h1 sol2_h2 sol3_h3

0.000 0.0000e+000 0.0000e+000 0.0000e+000 0.0000e+000 0.0000e+000 0.0000e+000 0.500 4.9659e-001 7.5910e-001 4.9659e-001 7.5910e-001 4.9659e-001 7.5910e-001 1.000 9.0431e-001 2.9880e+000 9.0430e-001 2.9882e+000 9.0430e-001 2.9882e+000

Risultati quesito 4: tabellaRisultati quesito 4: tabella

t solode45 (si riportano le 5 componenti)

0.000 0.0000e+000 1.0000e+000 0.0000e+000 0.0000e+000 5.0000e-001 0.000 2.0095e-004 1.0000e+000 -1.0106e-008 1.0056e-004 5.0080e-001 0.001 1.2057e-003 1.0000e+000 -3.6577e-007 6.0576e-004 5.0482e-001 0.006 6.2295e-003 1.0000e+000 -1.0024e-005 3.1924e-003 5.2492e-001 0.031 3.1348e-002 1.0000e+000 -2.8676e-004 1.7640e-002 6.2540e-001 0.131 1.3134e-001 9.9971e-001 -7.3360e-003 1.0019e-001 1.0258e+000 0.231 2.3124e-001 9.9801e-001 -2.9932e-002 2.2290e-001 1.4292e+000 0.331 3.3083e-001 9.9293e-001 -7.6340e-002 3.8634e-001 1.8419e+000 0.431 4.2963e-001 9.8165e-001 -1.5557e-001 5.9199e-001 2.2760e+000 0.531 5.2683e-001 9.6037e-001 -2.7846e-001 8.4293e-001 2.7523e+000 0.631 6.2120e-001 9.2402e-001 -4.5962e-001 1.1450e+000 3.3055e+000 0.731 7.1091e-001 8.6579e-001 -7.2085e-001 1.5085e+000 3.9948e+000 0.831 7.9331e-001 7.7604e-001 -1.0972e+000 1.9519e+000 4.9268e+000 0.931 8.6460e-001 6.4055e-001 -1.6488e+000 2.5087e+000 6.3116e+000 1.000 9.0430e-001 5.0996e-001 -2.1819e+000 2.9882e+000 7.7458e+000

t solode45 (si riportano le 5 componenti)

0.000 0.0000e+000 1.0000e+000 0.0000e+000 0.0000e+000 5.0000e-001 0.000 2.0095e-004 1.0000e+000 -1.0106e-008 1.0056e-004 5.0080e-001 0.001 1.2057e-003 1.0000e+000 -3.6577e-007 6.0576e-004 5.0482e-001 0.006 6.2295e-003 1.0000e+000 -1.0024e-005 3.1924e-003 5.2492e-001 0.031 3.1348e-002 1.0000e+000 -2.8676e-004 1.7640e-002 6.2540e-001 0.131 1.3134e-001 9.9971e-001 -7.3360e-003 1.0019e-001 1.0258e+000 0.231 2.3124e-001 9.9801e-001 -2.9932e-002 2.2290e-001 1.4292e+000 0.331 3.3083e-001 9.9293e-001 -7.6340e-002 3.8634e-001 1.8419e+000 0.431 4.2963e-001 9.8165e-001 -1.5557e-001 5.9199e-001 2.2760e+000 0.531 5.2683e-001 9.6037e-001 -2.7846e-001 8.4293e-001 2.7523e+000 0.631 6.2120e-001 9.2402e-001 -4.5962e-001 1.1450e+000 3.3055e+000 0.731 7.1091e-001 8.6579e-001 -7.2085e-001 1.5085e+000 3.9948e+000 0.831 7.9331e-001 7.7604e-001 -1.0972e+000 1.9519e+000 4.9268e+000 0.931 8.6460e-001 6.4055e-001 -1.6488e+000 2.5087e+000 6.3116e+000 1.000 9.0430e-001 5.0996e-001 -2.1819e+000 2.9882e+000 7.7458e+000

Risultati quesito 3: graficiRisultati quesito 3: grafici

Esercizio 3Esercizio 3Prova pratica del 26-06-06Prova pratica del 26-06-06

.141',171

],2,1()(')(2328

1)('' 3

yy

xxyxyxxy

Si consideri il seguente problema di Cauchy:

1- Si verifichi che la funzione è soluzione del problema proposto e si dica, motivando opportunamente la risposta, se tale soluzione è unica, sapendo che essa e la sua derivata sono funzioni limitate nell’intervallo [1, 2].

x

xxy162

Quesito 2Quesito 22- Si costruisca un file MATLABCognome_studente_matricola.m che, una volta avviato:a) faccia visualizzare una schermata con i dati personali ed una breve presentazione del problema;b) determini la soluzione approssimata utilizzando il metodo di Runge Kutta del quarto ordine considerando N1=100 e N2=200 intervalli; c) valuti l’errore assoluto nei nodi; d) faccia visualizzare una tabella riassuntiva in cui si riporti: intestazione: t soluzione1 soluzione2 errore1 errore2 e, ogni 5, i nodi ti , le soluzioni approssimate ed i

corrispondenti errori assoluti nei nodi coincidenti nei due casi, utilizzando i seguenti formati di stampa:

Quesiti 3 e 4Quesiti 3 e 4

3 cifre decimali e formato virgola fissa per i valori dei nodi,10 cifre decimali e virgola fissa per le soluzioni nei due casi,2 cifre decimali e formato floating point per gli errori.

3- Si esegua una figura con due finestre grafiche, una per N1=100 ed una per N2=200; in ciascuna finestra grafica si riporti la soluzione approssimata (color rosso e punto-linea) e la soluzione vera (color verde e linea continua).

4- Si commentino i risultati e si specifichi se essi soddisfano la aspettative teoriche.

Istruzioni quesito 2Istruzioni quesito 2clcclear all % Dati di input: t0=1;tmax=2;n1=100;n2=200;y0=[17 -14];f=strvcat('y(2)','1/8*(32+2*t.^3-y(1).*y(2))'); % Elaborazione:[T1,Y1]=Rungekutta4(t0,tmax,n1,y0,f);[T2,Y2]=RungeKutta4(t0,tmax,n2,y0,f);

xvera=T1.^2+16./T1;yvera=2*T1-16./T1.^2;vera=[xvera yvera];err1=abs(Y1-vera);err2=abs(Y2(1:2:end,:)-vera);

Istruzioni quesito 2Istruzioni quesito 2

tab=[T1 Y1 Y2(1:2:end,:) err1 err2];tab5=tab(1:5:end,:);tabella=[T1 Y1(:,1) Y2(1:2:end,1) err1(:,1) err2(:,1)];tabella5=tabella(1:5:end,:);

% Stampa della tabella:s='___________________________________________________';disp(s)fprintf(' t \t\t\t\t y1 y2 \t\t\t\t erry1 erry2 \n') disp(s);disp(' ')fprintf('%7.3f %14.10f %14.10f %12.2e %12.2e \n', tabella5')disp(' ')disp(s)

Istruzioni quesito 3Istruzioni quesito 3

% Stampa dei grafici:subplot(2,1,1) plot(T1,Y1(:,1),'.-r',T1,xvera,'g');grid onxlabel('tempo');ylabel('y1 e xvera');title(' RK4 (n1=100)') subplot(2,1,2)plot(T2,Y2(:,1),'.-r',T1,xvera,'g');grid onxlabel('tempo');ylabel('y2 e xvera');title(' RK4 (n2=200)')

Tabella risultatiTabella risultati__________________________________________________________________________________________________

t y1 y2 erry1 erry2 __________________________________________________________________________________ 1.000 17.0000000000 17.0000000000 0.00e+000 0.00e+000 1.050 16.3405952410 16.3405952383 2.95e-009 1.83e-010 1.100 15.7554545505 15.7554545458 5.07e-009 3.14e-010 1.150 15.2355434849 15.2355434787 6.61e-009 4.09e-010 1.200 14.7733333411 14.7733333338 7.73e-009 4.79e-010 1.250 14.3625000086 14.3625000005 8.56e-009 5.30e-010 1.300 13.9976923169 13.9976923083 9.16e-009 5.68e-010 1.350 13.6743518615 13.6743518524 9.62e-009 5.96e-010 1.400 13.3885714385 13.3885714292 9.95e-009 6.17e-010 1.450 13.1369827688 13.1369827593 1.02e-008 6.32e-010 1.500 12.9166666771 12.9166666673 1.04e-008 6.44e-010 1.550 12.7250806557 12.7250806458 1.05e-008 6.52e-010 1.600 12.5600000106 12.5600000007 1.06e-008 6.59e-010 1.650 12.4194697077 12.4194696976 1.07e-008 6.64e-010 1.700 12.3017647166 12.3017647065 1.08e-008 6.67e-010 1.750 12.2053571537 12.2053571435 1.08e-008 6.69e-010 1.800 12.1288888997 12.1288888896 1.08e-008 6.71e-010 1.850 12.0711486595 12.0711486493 1.08e-008 6.72e-010 1.900 12.0310526424 12.0310526323 1.09e-008 6.72e-010 1.950 12.0076282160 12.0076282058 1.09e-008 6.72e-010 2.000 12.0000000108 12.0000000007 1.08e-008 6.72e-010 ________________________________________________________________________________

GraficiGrafici

Si determini la traiettoria ed il tempo di attraversamento di un fiume largo 2 km, 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 4: Moto del BattelloEsercizio 4: Moto del Battello

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

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:

Eulero scalare:

Metodo di EuleroMetodo 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);

Soluzione del problema e simulazione Soluzione del problema e simulazione del motodel 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

L’estremo finale viene calcolato per tentativi.

Risultati con Risultati con tmax=0.65tmax=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

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 5: Problema Esercizio 5: Problema Preda-Preda-PredatorePredatore

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

Soluzione analitica del problemaSoluzione 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

Istruzioni: Metodo Eulero esplicito Istruzioni: Metodo Eulero esplicito e grafici di e grafici di PP( ( t t ),), QQ( ( t 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')

Grafici: Metodo di Eulero EsplicitoGrafici: Metodo di Eulero Esplicito

Istruzioni per grafico nel piano Istruzioni per grafico nel piano PP QQt0=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))])

Rappresentazione nel piano Rappresentazione nel piano PPQQ

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

1 0

2 0

d Pk P cPQ

d t

dQt k Q dPQ

d t

*

*

2 /

1/

P k d

Q k c

Istruzioni del Metodo HeunIstruzioni del Metodo Heun Grafici di Grafici di PP( ( t t ),), QQ( ( t 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')

Grafici prede, predatoriGrafici prede, predatoriMetodo di HeunMetodo di Heun

Soluzione nel piano Soluzione nel piano PPQQ::simulazione movimentosimulazione 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')

Simulazione: metodo di Heun

1 1

1 1

10

0 1 0 1 1

d xx t y t

dYd t t Y td td y

x t y td t

Yx y

Esercizio 6: Sistema differenziale Esercizio 6: Sistema differenziale lineare del 1° ordinelineare 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

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 6: quesiti c,d

a:a: esistenza, unicità della soluzioneesistenza, 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.

a: stabilità e condizionamentoa: 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.

n=input('n = '); % numero sottointervallit0=0;tmax=10;A=[-1 1; -1 -1];b=strvcat('0','0');y0=[-1 1];[T,Y]=Eulsis(t0,tmax,n,y0,A,b);str1='Eulero Implicito';

b: metodo di Eulero Implicitob: 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

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

% un problema differenziale di Cauchy di tipo vettoriale e le colonne

% di errore nella X e nella Y;

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

----------------------------------------------------------------------------------------------------- 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

tabella2.mtabella2.m : Eulero Implicito : Eulero Implicito n1 = 100n1 = 100

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

Errore per la Errore per la xx((tt): Eulero Implicito): Eulero Implicito

% Andamento degli errori nella valutazione di xplot(T,ErrX);grid title(’Andamento degli errori nella x')

tabella2.m:tabella2.m: Eulero Implicito Eulero Implicito n2 = 200n2 = 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

Soluzione con movimento nel piano Soluzione con movimento nel piano XYXY

n=200;t0=0;tmax=10;A=[-1 1; -1 -1];b=strvcat('0','0');y0=[-1 1];[T,Y]=Eulsis(t0,tmax,n,y0,A,b);plot(0,0,'or',y0(1),y0(2),'*g') % (0,0) punto % di stazionarietàhold onfor i=1:nplot(Y(i,1),Y(i,2),'ob')pause(0.25)title('Soluzione nel piano delle fasi')end

Grafico nel piano XYGrafico nel piano XY

Confronto risultati con n1=100 e n2=200Confronto risultati con n1=100 e n2=200clear allt0=0;tmax=10;y0=[-1 1];A=[-1 1; -1 -1];b=strvcat('0','0');n1=100;[T1,Y1]=Eulsis(t0,tmax,n1,y0,A,b);n2=200;[T2,Y2]=Eulsis(t0,tmax,n2,y0,A,b);t=T1; % calcolo degli erroriXv=exp(-t).*(sin(t)-cos(t));Yv=exp(-t).*(sin(t)+cos(t));ErrX1= abs(Y1(:,1)-Xv);ErrY1= abs(Y1(:,2)-Yv);% I vettori da inserire in tab devono essere della stessa

lunghezza e ErrX2,ErrY2 vanno calcolati negli stessi nodi di ErrX1, ErrY1!!!

ErrX2= abs(Y2(1:2:end,1)-Xv);ErrY2= abs(Y2(1:2:end,2)-Yv);

tab=[T1,ErrX1,ErrX2,ErrY1,ErrY2]; % costruzione della tabellas='---------------------------------------------------';disp(s)fprintf(' T ErrX1 ErrX2 ErrY1 ErrY2 \n');disp(s)fprintf('%6.3f %10.2e %10.2e %10.2e %10.2e\n',tab(1:10:end,:)')

Confronto dei risultatiConfronto 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-00510.000 4.85e-005 2.83e-005 3.47e-005 1.27e-005

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 4d: 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

Tabella Metodo RK4: Tabella Metodo RK4: n1 = 100n1 = 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

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

Errore per la x(t): Metodo RK4

Esercizio 7Esercizio 7Si consideri un satellite in orbita ellittica attorno ad un

pianeta di massa M. Supponiamo che il pianeta sia posto nell’origine del riferimento cartesiano bidimensionale e il problema sia normalizzato in modo che GM=1, con G costante di gravitazione universale.

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

322 4 aT

Quesiti 2a, 2b Quesiti 2a, 2b 1- Si 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 Runge-Kutta4 nel caso

In questo caso l’orbita è un’ellisse di centro (-1,0) e semiasse maggiore a = 2.

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

Quesiti 2c, 2d, 3Quesiti 2c, 2d, 3

c) faccia visualizzare una tabella riassuntiva, che riporti 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 2 finestre grafiche, si ritrovino per via grafica (corredata di label e titolo) le affermazioni del punto b) e si tracci l’andamento della soluzione in funzione del tempo.

3 - Si commentino i risultati.

Metodo RK4Metodo 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')

RisultatiRisultati

------------------------------------------------- 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

GraficiGrafici