1 / 48

7 a lezione - laboratorio

7 a lezione - laboratorio. Corso di Laurea ING. MECCANICA. a.a 2004-2005. Problema di Cauchy del 2° ordine nell’incognita. Esercizio 1: Problema del Pendolo. Si utilizzi il metodo di Runge Kutta 4 per calcolare la soluzione del problema del pendolo:.

kioshi
Download Presentation

7 a lezione - 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. 7a lezione - laboratorio Corso di Laurea ING.MECCANICA a.a 2004-2005

  2. Problema di Cauchy del 2° ordine nell’incognita Esercizio 1: Problema del Pendolo Si utilizzi il metodo di Runge Kutta 4 per calcolare la soluzione del problema del pendolo: g: accelerazione di gravità

  3. Problema equivalente del 1° ordine Sostituzioni Problema equivalente Risolvendo (*) in si ottengono: x (t)e quindi  (t), valore dell’angolo istante per istante, e , velocità.

  4. Dati del problema e metodo di soluzione Si assume: Si analizzano i risultati che si ottengono applicando il metodo di Runge Kutta 4.

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

  6. File “tabella.m” % Il file permette di scrivere una tabella in cui si riportano i valori % in uscita da un problema differenziale di Cauchy di tipo vettoriale. % Se si vuole riportare anche l'errore, vanno aggiunte altre colonne clc a=[T,x,y]; s='----------------------------------------------'; disp(s) fprintf('tempo\t\t\t angolo\t\t\t\t velocità\n') disp(s) n=length(T); for i=1:10:n fprintf('%6.3f %16.6e %16.6e \n',a(i,1:3)) end

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

  8. x  P y Movimento del punto P Il punto P ha coordinate: Per costruire la rappresentazione del punto P in movimento, determiniamo le coordinate di P, cioè: (sin ( x ( i )), - cos ( x ( i ))), i = 0,1,…,n. N.B. l’asse y è orientato verso il basso, cioè in modo contrario della normale rappresentazione.

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

  10. Esercizio 2: Moto del Battello Si determini la traiettoria ed il tempo di attraversamento di un fiume largo 2 km e con velocità della corrente di modulo s, da parte di un battello che si muove con velocità relativa (rispetto all’acqua) di modulo v, e che, partendo da un punto a valle (o a monte) del punto di attracco, si dirige sempre verso tale punto. Si utilizzi il metodo di Eulero. Si esaminerà il caso di partenza da un punto a valle.Lo studente può studiare l’altro caso.

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

  12. Modello del problema Velocità assoluta battello:

  13. [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:

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

  15. L’estremo finale viene calcolato per tentativi.

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

  17. Esercizio 3: 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.

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

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

  20. Grafici: Metodo di Eulero Esplicito

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

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

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

  24. Grafici prede, predatoriMetodo di Heun

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

  26. Simulazione: metodo di Heun

  27. Esercizio 4: 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.

  28. Esercizio 4: 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.

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

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

  31. b: metodo di Eulero Implicito n=input('n = '); % numero sottointervalli t0=0;tmax=10; coeff=[-1 1 0; -1 -1 0]; x0=-1;y0=1; [T,X,Y]=Eulsis(t0,tmax,n,x0,y0,coeff); str1='Eulero Implicito'; 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

  32. c: file “tabella2.m” % Il file scrive una tabella in cui si riportano i valori in uscita da % un problema differenziale di Cauchy di tipo vettoriale. % Si vuole riportare l'errore nella X e nella Y, quindi vanno % considerate le colonne relative. % ErrX= abs(X-Xv);ErrY= abs(Y-Yv); a=[T,X,Y, ErrX,ErrY]; s='------------------------------------------'; disp(s) fprintf(' T X Y ErrX ErrY \n'); disp(s) n=length(T); for i=1:10:n %stampa ogni 10 valori fprintf('%6.3f %16.6e %16.6e %10.2e %10.2e \n',a(i,:)) end

  33. Tabella Eulero Implicito: n1 = 100 ----------------------------------------------------------------------------------------------------- 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

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

  35. Errore per la x(t): Eulero Implicito

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

  37. Costruzione tabella riassuntiva n1=100; % numero sottointervalli t0=0;tmax=10;coeff=[-1 1 0; -1 -1 0]; x0=-1;y0=1; [T1,X1,Y1]=Eulsis(t0,tmax,n1,x0,y0,coeff); t=T1; Xv=exp(-t).*(sin(t)-cos(t)); Yv=exp(-t).*(sin(t)+cos(t)); ErrX1= abs(X1-Xv);ErrY1= abs(Y1-Yv); n2=200; % numero sottointervalli [T2,X2,Y2]=Eulsis(t0,tmax,n2,x0,y0,coeff); t=T2; Xv=exp(-t).*(sin(t)-cos(t)); Yv=exp(-t).*(sin(t)+cos(t)); ErrX2= abs(X2-Xv);ErrY2= abs(Y2-Yv); a=[T1,ErrX1,ErrX2(1:2:end),ErrY1,ErrY2(1:2:end)]; fprintf(' T ErrX1 ErrX2 ErrY1 ErrY2\n'); fprintf('%6.3f %10.2e %10.2e %10.2e %10.2e \n',a(1:10:end,:)')

  38. Tabella riassuntiva dei risultati T ErrX1 ErrX2 ErrY1 ErrY2 0.000 0.00e+000 0.00e+000 0.00e+000 0.00e+000 1.000 4.75e-002 2.46e-002 1.11e-002 5.55e-003 2.000 1.36e-002 6.75e-003 3.33e-002 1.73e-002 3.000 1.06e-002 5.80e-003 1.68e-002 8.42e-003 4.000 9.74e-003 5.02e-003 5.17e-004 4.12e-005 5.000 2.79e-003 1.28e-003 3.51e-003 1.92e-003 6.000 5.50e-004 3.86e-004 1.90e-003 9.42e-004 7.000 7.99e-004 4.28e-004 2.79e-004 8.94e-005 8.000 2.90e-004 1.32e-004 2.05e-004 1.28e-004 9.000 2.80e-007 1.32e-005 1.46e-004 7.47e-005 10.000 4.85e-005 2.83e-005 3.47e-005 1.27e-005

  39. MODIFICHE ai files EULSIS e TRAPEZI .. …. .. . %file TRAPEZI.M while t<tmax c1_old=eval(coeff(1,:)); c2_old=eval(coeff(2,:)); t=t+h; if abs(tmax-t)<1.e-13 t=tmax; end c1_new=eval(coeff(1,:)); c2_new=eval(coeff(2,:)); .. …. .. . end .. …. .. . %file EULSIS.M while t < tmax tn=x1+h*tn1; x1=mat\tn; X=[X;x1(1)]; Y=[Y;x1(2)]; t=t+h; if abs(tmax-t)<1.e-13 t=tmax; end T=[T;t]; end

  40. Soluzione con movimento nel piano XY n=200; t0=0; tmax=10; coeff=[-1 1 0;-1 -1 0]; x0=-1;y0=1; [T,X,Y]=Eulsis(t0,tmax,n,x0,y0,coeff); plot(0,0,'or',x0,y0,'*g') % (0,0) punto % di stazionarietà hold on for i=1:n plot(X(i),Y(i),'ob') pause(0.25) end

  41. Grafico nel piano XY

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

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

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

  45. Errore per la x(t): Metodo RK4

More Related