1 / 68

11 a lezione di laboratorio

11 a 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:.

kimama
Download Presentation

11 a lezione di laboratorio

An Image/Link below is provided (as is) to download presentation Download Policy: Content on the Website is provided to you AS IS for your information and personal use and may not be sold / licensed / shared on other websites without getting consent from its author. Content is provided to you AS IS for your information and personal use only. Download presentation by click this link. While downloading, if for some reason you are not able to download a presentation, the publisher may have deleted the file from their server. During download, if you can't get a presentation, the file might be deleted by the publisher.

E N D

Presentation Transcript


  1. 11a lezione di laboratorio Laurea Ingegneria CIVILE Lauree Specialistiche in Ingegneria CHIMICA, ELETTRONICA, AMBIENTE a.a. 2007-2008

  2. Esercizio 1Prova pratica del 04-09-2007 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].

  3. 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; • faccia visualizzare una tabella riassuntiva in cui si • riporti:

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

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

  6. 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)';

  7. 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')

  8. Risultati 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

  9. Risultati quesito 3: grafici

  10. Esercizio 2Prova pratica del 30-11-2006 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.

  11. Quesito 2 2- Si costruisca un file MATLAB: Cognome_studente_matricola.mche una volta avviato: • 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.

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

  13. Riduzione a sistema del primo ordine Si ottiene quindi il sistema:

  14. 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,:);

  15. 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)')

  16. Risultati quesito 2: tabella tsol1_h1sol2_h2sol3_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

  17. Risultati 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

  18. Risultati quesito 3: grafici

  19. Esercizio 3Prova pratica del 26-06-06 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].

  20. Quesito 2 • 2- Si costruisca un file MATLAB • Cognome_studente_matricola.m che, una volta avviato: • 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:

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

  22. 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);

  23. 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)

  24. 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)')

  25. 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 ________________________________________________________________________________

  26. Grafici

  27. 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 puntoa valle.Lo studente può studiare l’altro caso.

  28. 2 km s s P P y y V V dy/dt dy/dt dx/dt dx/dt 0.1 km v v x x A A

  29. Modello del problema Velocità assoluta battello:

  30. [T,Y]=Eulero(t0,tmax,n,y0,f) f=strvcat(f1,f2); [T,Y]=Eulero(t0,tmax,n,y0,f) [T,X,Y]=Eulero2(t0,tmax,n,x0,y0,f1,f2); Metodo di Eulero Eulero scalare: Eulero vettoriale: Eulero vettoriale applicato in serie:

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

  32. L’estremo finale viene calcolato per tentativi.

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

  34. Esercizio 5: Problema Preda-Predatore Sistema differenziale non lineare di ordine 1 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.

  35. Soluzione analitica del problema Variabili separate Integrando i due membri si ottiene l’Int. Gen. cost. si calcola imponendo le condizioni iniziali.

  36. 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')

  37. Grafici: Metodo di Eulero Esplicito

  38. Istruzioni per grafico nel piano PQ 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))])

  39. Rappresentazione nel piano PQ Il punto segnato con * è il punto critico che si ottiene ponendo:

  40. 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')

  41. Grafici prede, predatoriMetodo di Heun

  42. 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')

  43. Simulazione: metodo di Heun

  44. Esercizio 6: Sistema differenziale lineare del 1° ordine • Stabilire se il problema ammette soluzione unica e se è ben condizionato. • Calcolare la soluzione con il metodo di Eulero Implicito,n1=100 e n2=200; • n1, n2= numero sottointervalli.

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

  46. a:esistenza, unicità della soluzione Esiste unica la soluzione del problema.

  47. a: stabilità e condizionamento La matrice dei coefficienti è i suoi autovalori si trovano risolvendo l’equazione caratteristica: il sistema è asintoticamente stabile e quindi ben condizionato.

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

  49. 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,:)')

More Related