a
11
lezione di laboratorio
Laurea Ingegneria CIVILE
Lauree Specialistiche in Ingegneria
CHIMICA, ELETTRONICA,
AMBIENTE
a.a. 2007-2008
Esercizio 1
Prova pratica del 04-09-2007
Si consideri il seguente problema di Cauchy:


1  y( x ) y' ( x ) y' ' ( x )  1  y'2 ( x ) y' ( x )

y' 0  2.
 y0  2,
x  (0,2],
1- Si verifichi che la funzione y x   2e x è 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].
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) 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 errRK
in 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 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:
variabili di input delle functions
clear all
clc
t0=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:
applicazione 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 grafici
fprintf(' 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: tabella
nodi
0.000
0.200
0.400
0.600
0.800
1.000
1.200
1.400
1.600
1.800
2.000
solH
2.0000000000
2.4426790826
2.9833405503
3.6436717793
4.4501604196
5.4351568857
6.6381720177
8.1074619672
9.9019638802
12.0936600235
14.7704651859
solRK
2.0000000000
2.4428055062
2.9836493707
3.6442375557
4.4510817835
5.4365635447
6.6402336811
8.1103996994
9.9060645217
12.0992944794
14.7781115880
errH
0.00e+000
5.18e-005
1.04e-004
1.55e-004
2.07e-004
2.59e-004
3.11e-004
3.62e-004
4.14e-004
4.66e-004
5.17e-004
errRK
0.00e+000
4.13e-009
8.25e-009
1.24e-008
1.65e-008
2.06e-008
2.48e-008
2.89e-008
3.30e-008
3.71e-008
4.13e-008
Risultati quesito 3: grafici
Esercizio 2
Prova pratica del 30-11-2006
Si consideri il problema differenziale di Cauchy
 y' ' '  x 2 y y' ' y z' , x  (0,1]

 z' '  x z z'4 y'
 y(0)  0, y' (0)  1, y' ' (0)  0,

1
 z (0)  0, z' (0)  .
2

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
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 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.
Riduzione a sistema del primo ordine
y  y  1 , y   1  y  2  , y   2   y  3 
z  y  4  , y  4   y  5 
Si ottiene quindi il sistema:

y  1  y  2 


y
 2  y  3


2

y
3

x
y  1 y  3   y  1 y  5 




y  4   y  5 

 y  5   x y  4  y  5   4 y  2 
Risoluzione quesito 2
clear all;clc
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])];
% prima e quarta colonna della soluzione numerica
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'];
end
tab=[T(1:step:end) tab];tabr=tab(1:5:end,:);
Risoluzione quesiti 3 e 4
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)];
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)')
Risultati quesito 2: tabella
t
0.000
0.500
1.000
sol1_h1
0.0000e+000
4.9659e-001
9.0431e-001
0.0000e+000
7.5910e-001
2.9880e+000
sol2_h2
0.0000e+000
4.9659e-001
9.0430e-001
0.0000e+000
7.5910e-001
2.9882e+000
sol3_h3
0.0000e+000
4.9659e-001
9.0430e-001
0.0000e+000
7.5910e-001
2.9882e+000
Risultati quesito 4: tabella
t
0.000
0.000
0.001
0.006
0.031
0.131
0.231
0.331
0.431
0.531
0.631
0.731
0.831
0.931
1.000
solode45 (si riportano le 5 componenti)
0.0000e+000
2.0095e-004
1.2057e-003
6.2295e-003
3.1348e-002
1.3134e-001
2.3124e-001
3.3083e-001
4.2963e-001
5.2683e-001
6.2120e-001
7.1091e-001
7.9331e-001
8.6460e-001
9.0430e-001
1.0000e+000
1.0000e+000
1.0000e+000
1.0000e+000
1.0000e+000
9.9971e-001
9.9801e-001
9.9293e-001
9.8165e-001
9.6037e-001
9.2402e-001
8.6579e-001
7.7604e-001
6.4055e-001
5.0996e-001
0.0000e+000
-1.0106e-008
-3.6577e-007
-1.0024e-005
-2.8676e-004
-7.3360e-003
-2.9932e-002
-7.6340e-002
-1.5557e-001
-2.7846e-001
-4.5962e-001
-7.2085e-001
-1.0972e+000
-1.6488e+000
-2.1819e+000
0.0000e+000
1.0056e-004
6.0576e-004
3.1924e-003
1.7640e-002
1.0019e-001
2.2290e-001
3.8634e-001
5.9199e-001
8.4293e-001
1.1450e+000
1.5085e+000
1.9519e+000
2.5087e+000
2.9882e+000
5.0000e-001
5.0080e-001
5.0482e-001
5.2492e-001
6.2540e-001
1.0258e+000
1.4292e+000
1.8419e+000
2.2760e+000
2.7523e+000
3.3055e+000
3.9948e+000
4.9268e+000
6.3116e+000
7.7458e+000
Risultati quesito 3: grafici
Esercizio 3
Prova pratica del 26-06-06
Si consideri il seguente problema di Cauchy:

1

 y' ' ( x )  32  2 x 3  y( x ) y' ( x )

8
 y1  17,
y' 1  14.

x  (1,2],
16
y x   x 
1- Si verifichi che la funzione
x è 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].
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) 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 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 2
clc
clear 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 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 3
% Stampa dei grafici:
subplot(2,1,1)
plot(T1,Y1(:,1),'.-r',T1,xvera,'g');grid on
xlabel('tempo');ylabel('y1 e xvera');
title(' RK4 (n1=100)')
subplot(2,1,2)
plot(T2,Y2(:,1),'.-r',T1,xvera,'g');grid on
xlabel('tempo');ylabel('y2 e xvera');
title(' RK4 (n2=200)')
Tabella 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
________________________________________________________________________________
Grafici
Esercizio 4: Moto del Battello
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.
2 km
s
y
P
V
dy/dt
0.1 km
dx/dt
v
x
A
Modello del problema
 dx


 dt

 dy
 
 dt

x 0   2
vx
x y
vy
2
x y
2
2
2
s
Velocità assoluta battello:
y 0   0.1
dx
dy
V 
i
j
dt
dt
Metodo di Eulero
Eulero scalare:
z i 1  z i  h * f t i , z i 
[T,Y]=Eulero(t0,tmax,n,y0,f)
Eulero vettoriale:
 x i 1  x i  h * f 1 t i , x i , y i 

 y i 1  y i  h * f 2 t i , x i , y i 
f=strvcat(f1,f2);
[T,Y]=Eulero(t0,tmax,n,y0,f)
Eulero vettoriale applicato in serie:
 x i 1

 y i 1
 x i 1  x i  h * y i
 x i  h * f 1 t i , x i , y i 

  y i 1  y i  h * g sin x i 1 
 y i  h * f 2 t i , x i 1 , y i  
 x 0  0; y 0  0.5
[T,X,Y]=Eulero2(t0,tmax,n,x0,y0,f1,f2);
Soluzione del problema e simulazione
del moto
t0=0;tmax=0.5; %primo valore di tentativo
y0=[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 on
for i=1:length(T)
plot(x(i),y(i),'*b') %*=simbolo,b=blu
pause(0.30)
end
hold off
title([' Traiettoria del battello - tmax = ' num2str(tmax)])
xlabel(' X [km]');ylabel(' Y [km]')
grid
L’estremo finale viene calcolato per tentativi.
Risultati con tmax=0.65
T
X
0.000
2.000000e+000
0.013
1.935081e+000
0.026
1.870240e+000
0.039
1.805508e+000
0.052
1.740916e+000
... ... ...
0.559
1.919500e-002
0.572
1.238269e-002
0.585
7.255724e-003
0.598
3.652893e-003
0.611
1.388678e-003
0.624
2.439726e-004
0.637
-5.607722e-005
0.650
7.966688e-005
Y
1.000000e-001
1.357541e-001
1.702052e-001
2.033141e-001
2.350406e-001
1.821415e-001
1.564994e-001
1.307019e-001
1.048019e-001
7.884131e-002
5.285139e-002
2.685208e-002
8.522218e-004
Esercizio 5: Problema
Preda-Predatore
 P t 

 Q t 
 P  0  P
0

k1 P  t   c P  t  Q  t 
 k2Q t   d P t Q t 
Q  0   Q0
Sistema differenziale non lineare di ordine 1
P(t): prede, Q(t): predatori.
Si trovi la soluzione del problema per t  0, 3
assumendo: k1 = 2; k2 = 10; c = 0.001;
d = 0.002; P0=5000; Q0=100.
Soluzione analitica del problema
dQ k 2Q  d PQ  k 2  d P  Q  k 2  d P 
Q



dP
k1P  c PQ
P
 k1  c Q  P
 k1  c Q 
k1  c Q
k 2  d P
dQ 
dP Variabili separate
Q
P
Integrando i due membri si ottiene l’Int. Gen.
k1ln Q  c Q   k 2 ln P  d P  cost.
 k1ln Q  c Q  k 2 ln P  d P  cost.
cost. si calcola imponendo le condizioni
iniziali.
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’); grid
title(['Risultati del metodo di ',str1])
legend('prede','predatori')
title(['Risultati del metodo di ',str1,...
' con h =' num2str(h)])
xlabel('Tempo')
Grafici: Metodo di Eulero Esplicito
Istruzioni per grafico nel piano P Q
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)+ .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 PQ
Il punto segnato
con * è il punto
critico che si
ottiene ponendo:
 dP
*
 d t  k 1P  cPQ  0

P  k2 / d

 *

 dQ  t    k 2Q  dPQ  0
 Q  k1/ c
 d t
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)); grid
legend('prede','predatori')
title(['Risultati del metodo di ',str1,...
' con h = ',num2str(h)])
xlabel('Tempo')
Grafici prede, predatori
Metodo di Heun
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')
Simulazione: metodo di Heun
Esercizio 6: Sistema differenziale
lineare del 1° ordine
 dx
 x t   y t 
 dt 
 dY
 1 1 
Y  t 

t   


 dy
 dt
 1 1

 x t   y t   

  1
 dt

Y  0   
 x  0   1

y  0  1
1




t
0,10
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.
Esercizio 6: quesiti c,d
c) Si calcoli l’errore nei nodi
soluzione vera è:
t i sapendo che la
 x  t   e  sint  cost 

t
 y  t   e  sint  cost 
t
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.
a: esistenza, unicità della soluzione

 x t    x t   y t   f


 y t    x t   y t   f
 t , x  t  , y  t    Y t  F t ,Y
  
 t , x  t  , y  t 
1
2
 Y
 
F
J
 1 1 

,

Y
  1  1
F
J
2


Esiste unica la soluzione del problema.
a: stabilità e condizionamento
 1 1 
A


1

1


La matrice dei coefficienti è
i
suoi autovalori si trovano risolvendo l’equazione
caratteristica:
 1    1     1  0  1,2  1
i  e  1,2   0

il sistema è asintoticamente stabile e quindi
ben condizionato.
b: metodo di Eulero Implicito
n=input('n = '); % numero sottointervalli
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);
str1='Eulero Implicito';
Stampa dei risultati
disp(['Risultati del metodo di ',str1]);
t=T;
Xv=exp(-t).*(sin(t)-cos(t));
Yv=exp(-t).*(sin(t)+cos(t));
tabella2
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 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 valori
fprintf('%6.3f %16.6e %16.6e %10.2e %10.2e\n',
tabella(1:10:end,:)')
tabella2.m : Eulero Implicito
n1 = 100
----------------------------------------------------------------------------------------------------T
X
Y
ErrX
ErrY
----------------------------------------------------------------------------------------------------0.000
1.000
2.000
3.000
4.000
5.000
6.000
7.000
8.000
9.000
10.000
-1.000000e+000
6.326408e-002
1.657577e-001
6.695044e-002
7.847666e-003
-5.585724e-003
-3.622299e-003
-8.876485e-004
9.098508e-005
1.630219e-004
6.190752e-005
1.000000e+000
5.194194e-001
1.000365e-001
-2.547604e-002
-2.531598e-002
-8.060362e-003
-2.110340e-004
1.007195e-003
4.883277e-004
8.486890e-005
-2.813839e-005
0.00e+000
4.75e-002
1.36e-002
1.06e-002
9.74e-003
2.79e-003
5.50e-004
7.99e-004
2.90e-004
2.80e-007
4.85e-005
0.00e+000
1.11e-002
3.33e-002
1.68e-002
5.17e-004
3.51e-003
1.90e-003
2.79e-004
2.05e-004
1.46e-004
3.47e-005
% Grafico componente x della soluzione vera ed
approssimata
plot(T,Xv,T,Y(:,1));grid
xlabel('Tempo (s)')
title(’Metodo di Eulero implicito - n=100')
legend('Xvera','Xapp')
Errore per la x(t): Eulero Implicito
% Andamento degli errori nella valutazione di x
plot(T,ErrX);grid
title(’Andamento degli errori nella x')
tabella2.m: 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
Soluzione con movimento nel piano
XY
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 on
for i=1:n
plot(Y(i,1),Y(i,2),'ob')
pause(0.25)
title('Soluzione nel piano delle fasi')
end
Grafico nel piano XY
Confronto risultati con n1=100 e n2=200
clear all
t0=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 errori
Xv=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 tabella
s='---------------------------------------------------';
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 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
d: metodo di Runge-Kutta 4
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';
Stampa dei risultati
disp(['Risultati del metodo di ',str1]);
t=T;
Xv=exp(-t).*(sin(t)-cos(t));
Yv=exp(-t).*(sin(t)+cos(t));
tabella
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
% Grafico soluzione vera ed approssimata
plot(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 7
Si 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:
 d 2x
2
2 3/ 2


x
/(
x

y
)
 2
 dt
 2
 d y   y /( x 2  y 2 )3 / 2 .
 d t 2
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 T 2  4 2 a 3
dove a è il semiasse maggiore dell’orbita.
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
x(0)  1,
y(0)  0,
x '(0)  0,
y '(0)  6 / 2;
In questo caso l’orbita è un’ellisse di centro (-1,0)
e semiasse maggiore a = 2.
Quesiti 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 RK4
t0=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')
Risultati
------------------------------------------------tempo
soluzione
derivata_soluzione
------------------------------------------------0.000
1.777
3.554
5.331
7.109
8.886
10.663
12.440
14.217
15.994
17.772
1.00000000
- 0.03265426
- 1.35404660
- 2.28448065
- 2.82363255
- 2.99985645
- 2.82341431
- 2.28404148
- 1.35338881
- 0.03186270
1.00000057
0.00000000
1.51593990
1.70460966
1.32746219
0.71084577
-0.00027980
-0.71137187
-1.32787730
-1.70479154
-1.51561630
0.00109146
0.00000000
-0.81632406
-0.63935333
-0.41023654
-0.19934779
0.00006125
0.19947120
0.41036158
0.63947005
0.81630315
-0.00092295
1.22474487
0.39064849
-0.09962284
-0.29773434
-0.38356134
-0.40826680
-0.38352275
-0.29764433
-0.09943831
0.39107078
1.22473755
Grafici