7 a- 8 a lezione di laboratorio

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

2. Comando ezsurf figure(1) z='X.*exp(-(X.^2+Y.^2))'; ezsurf(z,[-2,2,-2,2]); %se non si specifica l’insieme la superficie è disegnata nel dominio di default -2*pi<x<2*pi,-2*pi<y<2*pi colorbar; title(' ezsurf')

3. Comando ezcontour Con il comando contour tracciamo le linee di livello nel dominio fissato, se non si fornisce vengono plottate nel dominio di default figure(2) z='X.*exp(-(X.^2+Y.^2))'; ezcontour(z,[-2,2,-2,2]);

4. Comando ezsurf Il comando ezsurfpermette di rappresentare anche superfici date in coordinate parametriche ad esempio: figure(3) funx='2*cos(s)'; funy='2*sin(s)'; funz='z'; ezsurf(funx,funy,funz)

5. Comando surf % Le istruzioni servono per i tre grafici che seguono. x=-2:.2:2;y=-2:.2:2; [X,Y]=meshgrid(x,y); Z=X.*exp(-(X.^2+Y.^2)); % comando surf figure(4) surf(X,Y,Z);colorbar title('surf')

6. Comando contour figure(5) contour(X,Y,Z,20) % si specifica il numero di curve %contour(X,Y,Z,[-.4:.2:.4]) %si specificano i valori in cui si vogliono le curve title('contour')

7. Comando quiver figure(4) [px,py]=gradient(Z,.2,.2);%[px,py]=gradient(Z); quiver(X,Y,px,py) title('quiver')

8. quiver e contour figure(5) contour(X,Y,Z,20;hold on quiver(X,Y,px,py);hold off

9. Esercizio 1 Sia dato il seguente problema alle derivate parziali (pde):

10. Quesiti a, b a) Si verifichi che la funzione: è soluzione del problema. • Si valuti l’errore assoluto che si commette se si usa il “metodo upwind” ed il “metodo implicito”, fissando il numero di intervalli temporali M = 10, al variare del passo temporale k e considerando il valore del passo spaziale h=0.25. • Si indichi con N il numero degli intervalli spaziali sull’asse x.

11. Soluzione del quesito a): Verifica Quindi e:

12. Quesito b): Metodo UPWIND Approssimazioni utilizzate: Indicando quindi si ottiene:

13. Problema discreto il metodo CONVERGE quando . Se si assume per ogni livello temporale j: il problema discreto diventa: Sappiamo che se:

14. Costruzione delle formule Dalle relazioni: si ottiene, per ogni livello temporale j, tenendo anche conto della condizione al contorno la seguente equazione vettoriale:

15. Forma della matrice A Si deduce allora che la matrice A, di dimensioniNxN ,per ogni j, ha la forma: e, nell’ipotesi risulta: Il metodo è condizionatamente stabile!!!

16. Quesito b): Metodo IMPLICITO In questo caso si colloca la pde in e si approssimano le derivate parziali con: due differenze all’indietro!! Si indica ancora: il problema diventa:

17. Problema discretizzato con il metodo implicito Indicando ancora: si ottiene il sistema:

18. Sistema relativo al metodo implicito il metodo converge quindi per: La formula per l’errore è: Poiché risulta: Convergenza incondizionata!!!

19. Istruzioni relative al quesito b) clear all; clc t0=0;x0=0;xN=10;h=0.25; M=10;c='3';c1=eval(c); f='(x-2).*exp(-2*(x-2).^2)';% cond. iniziale g='-(3*t+2).*exp(-2*(3*t+2).^2)';%cond.contorno r='0';%termine noto fprintf(['M =',num2str(M),'\n\n h k k+h alpha err_imp err_up \n']) Uveras='(X-c1*T-2).*exp(-2*(X-c1*T-2).^2)'; for k=[0.05 h/3 0.1 0.5] alpha=c1*k/h; hpk=h+k; [x,t,sol1]=PDE_upwind(t0,M,x0,xN,h,k,c,r,f,g); [x,t,sol2]=PDE_implicito(t0,M,x0,xN,h,k,c,r,f,g); [X,T] = meshgrid(x,t);Uvera=eval(Uveras); err1=abs(Uvera-sol1);% matrice degli errori: upwind err2=abs(Uvera-sol2);% matrice degli errori: implicito errore_max_up=max(max(err1)); errore_max_imp=max(max(err2)); tab=[h k h+k alpha errore_max_imp errore_max_up]; fprintf('%6.2f %8.4f %8.4f %6.2f %13.4e %13.4e \n',tab') end

20. Function PDE_upwind x=(x0:h:xN)'; x(end)=xN; N=length(x)-1; tM=k*M+t0; t=linspace(t0,tM,M+1)'; U0=eval(f).*ones(N+1,1); %condizione iniziale U(x,t0) vv=eval(g).*ones(M+1,1); %condizione al contorno U(x0,t) Vj=zeros(N,1); Uj=U0(2:N+1);sol=U0'; t=t0;x=x(2:end); for j=1:M alpha=(eval(c)*k/h).*ones(N,1); tnoto=eval(r).*ones(N,1); A=diag(1-alpha)+diag(alpha(2:end),-1); Vj(1)=vv(j); Uj1=A*Uj+alpha(1)*Vj +k*tnoto; sol=[sol;[vv(j+1); Uj1]']; Uj=Uj1; t=t+k; end t=linspace(t0,tM,M+1)';x=[x0;x];

21. Function PDE_implicito x=(x0:h:xN)'; x(end)=xN; N=length(x)-1; tM=k*M+t0; t=linspace(t0,tM,M+1)'; U0=eval(f).*ones(N+1,1); %condizione iniziale U(x,t0) vv=eval(g).*ones(M+1,1);%condizione al contorno U(x0,t) Vj1=zeros(N,1);sol=U0';Uj=U0(2:N+1); t=t0;x=x(2:end); for j=1:M t=t+k; alpha=(eval(c)*k/h).*ones(N,1); tnoto=eval(r).*ones(N,1); A=-diag(alpha(2:end),-1)+diag(1+alpha); Vj1(1)=vv(j+1); b=Uj+alpha(1)*Vj1 +k*tnoto; Uj1=A\b; sol=[sol;[vv(j+1); Uj1]']; Uj=Uj1; end t=linspace(t0,tM,M+1)';x=[x0;x];

22. Risultati al variare del passo k entrambi i metodi sono consistenti: sono stabili, e quindi convergono. l’ implicito converge, upwind è instabile! e quindi non converge. M = 10

23. Osservazione sul caso Linee caratteristiche : il metodo upwind fornisce:

24. Commenti sul caso Si ottiene, per la forma di Sono valori corretti perché assegnati dalle condizioni. Lo stesso per j > 0. In questo caso, il metodo upwind calcola la soluzione esatta, i nodi sono tutti sulle rette caratteristiche!!!!

25. Rappresentazione della soluzione e delle curve di livello % % Rappresentazione della superficie e delle % curve di livello % k=h/3 [X,T]=meshgrid(x,t); figure(1) S=surfl(X,T,sol1); %surfl title('soluzione approssimata:metodo upwind') xlabel('x'),ylabel('t') figure(2) C=contour(X,T,sol1,20); %20 curve di livello title('Curve di livello') xlabel('x'),ylabel('t')

26. Superficie: metodo upwind h=0.25, k=h/3

27. Andamento della soluzione al variare di t per x fissato. Si ottiene selezionando Figure Palette dal menu del tasto View; sulla sinistra compare la lista delle variabili coinvolte. La figura presentata si ottiene premendo su sol1. Cliccando su una linea si individua a quale componente della sol1 corrisponde.

28. Migliore definizione dei comandi PLOT, SURF, CONTOUR Se si vuole definire meglio le figure, conviene utilizzare istruzioni del tipo: H=surf(X,T,sol1); set(gca,'Fontsize',14) % 14 punti per pollice set(H, 'LineWidth',2) % spessore della linea Istruzioni analoghe per plot e contour

29. Esercizio 2 Si consideri il seguente problema misto ai valori iniziali ed al contorno, con coefficienti non costanti:

30. Quesiti a, b • Si determinino le linee caratteristiche e si • verifichi che la soluzione del problema ai • valori iniziali su tutto l’asse reale, soddisfa anche la condizione al bordo per x = 0. b) Si valuti il massimo dell’errore assoluto che si commette usando il “metodo upwind” ed il “metodo implicito” se si fissa il tempo finale tM=3 e si prendono i passi spaziali h=0.5,0.2,0.1.

31. Soluzione del quesito a) Per individuare le caratteristiche della pde data, si risolve il problema di Cauchy: Esso ha la soluzione e p(x,t) è la linea che collega Si verifichi ora che è la soluzione del problema pde+ condizione iniziale che soddisfa anche la condizione assegnata al bordo.

32. Soluzione del quesito b) clear all; clc t0=0;tM=3; % in questo caso si assegna tmax x0=0;xN=3;c='2*t';t=tM;c1=eval(c);h=[0.2 0.1 0.05]'; k=h./c1; M=round((tM-t0)./k); alpha=c1*k./h; f='1./(1+x.^2)'; %condizione iniziale g='1./(1+t.^4)'; %condizione al contorno r='0'; Uveras='1./(1+(X-T.^2).^2)'; % soluzione vera tab=[]; for i=1:length(h) [x,t,sol1]=PDE_upwind(t0,M,x0,xN,h(i),k(i),c,r,f,g); [x,t,sol2]=PDE_implicito(t0,M,x0,xN,h(i),k(i),c,r,f,g); %soluzione vera e errore massimo del metodo [X T]=meshgrid(x,t);Uvera=eval(Uveras); if i==1 %grafici per h=0.2 e k=h/6 grafici end err1=abs(Uvera-sol1);err2=abs(Uvera-sol2); err1max=max(max(err1));err2max=max(max(err2)); tab=[tab;err2max err1max]; end tab=[h k h+k alpha tab]; fprintf(['h k k+h alpha err_imp err_upw \n']) fprintf('%6.2f %8.4f %8.4f %6.2f %13.4e %13.4e \n',tab')

33. File grafici (prima parte) figure() surf(X,T,Uvera) set(gca, 'FontWeight','bold','Fontsize',12) title('Soluzione vera');xlabel('x');ylabel('t') titolo1=['- h =', num2str(h(i))]; for m=1:2 if m==1 sol=sol1; titolo=['metodo upwind',titolo1]; elseif m==2 sol=sol2; titolo=['metodo implicito',titolo1]; end

34. File grafici (seconda parte) figure() surf(X,T,sol) set(gca, 'FontWeight','bold','Fontsize',12) title(['Soluzione approssimata:', titolo]); xlabel('x');ylabel('t') figure() [C,H]=contour(X,T,sol,20);% 20 linee di livello set(gca, 'FontWeight','bold','Fontsize',12) set(H,'LineWidth',2) title(['Curve di livello:',titolo]) xlabel('x'); ylabel('t')

35. Rappresentazione della soluzione vera

36. Superficie approssimata: metodo upwind k=h/6

37. Curve di livello: metodo upwind k=h/6

38. Superficie approssimata: metodo implicito k=h/6

39. Curve di livello: metodo implicito k=h/6

40. Errori in t= tM=3 h k k+h alpha err_imp err_upw 0.20 0.0333 0.2333 1.00 3.1685e-001 2.0392e-001 0.10 0.0167 0.1167 1.00 2.1497e-001 1.2541e-001 0.05 0.0083 0.0583 1.00 1.3615e-001 7.2635e-002 La tabella si riferisce al tempo finale tM=3; i valori di k sono stati calcolati con la relazione k=h/c(tM) dove c(tM)=2*tM e quindik=h/6.

41. Esercizio 3 Sia dato il seguente problema alle derivate parziali a coefficienti non costanti: con soluzione vera:

42. Quesiti 1) e 2) • Si verifichi che la funzione (1) è soluzione • del problema proposto e si calcoli, in • corrispondenza del passo spaziale h=0.2, • il passo temporale k massimo per cui il • metodo esplicito converge. 2) Si valuti, per il passo spaziale h=0.2 e fissando il tempo finale tM=1, l’errore assoluto massimo che si commette usando il “metodo upwind”ed il “metodo implicito”.

43. Quesito 3) 3) Si costruisca una tabella che riporti l’intestazione: t sol1 sol2 err1 err2 con le quantità t, sol1, sol2, err1, err2 rappresentanti rispettivamente, i nodi temporali, la soluzione numerica e l’errore ottenuti con i due metodi, da riportare uno ogni due, valutati in corrispondenza del valore x=2,utilizzando i seguenti formati di stampa: 3 cifre decimali e formato virgola fissa per i nodi, 6 cifre decimali e formato esponenziale per la soluzione nei due metodi, 2 cifre decimali e formato virgola mobile per l’errore nei due metodi.

44. Istruzioni per risolvere i quesiti 1) e 2) clear all; clc t0=0;tM=1; x0=0;xN=4; h=0.2; c='t.^2+x'; t=tM;x=xN;k=h/eval(c);M=round((tM-t0)/k); r='exp(-t).*(2*t.^2+x).*x'; f='x.^2'; % condizione iniziale U(x,t0) g='0'; % condizione al contorno U(x0,t) [x,t,sol1]=PDE_upwind(t0,M,x0,xN,h,k,c,r,f,g); [x,t,sol2]=PDE_implicito(t0,M,x0,xN,h,k,c,r,f,g); [X T]=meshgrid(x,t); Uvera=X.^2.*exp(-T); % soluzione vera err1=abs(Uvera-sol1); err2=abs(Uvera-sol2);

45. Costruzione delle tabelle: quesiti 2) e 3) err1max=max(max(err1));% massimo dell’errore err2max=max(max(err2));% massimo dell’errore tab=[h k h+k err2max err1max]; fprintf([h k k+h err_imp err_upw \n']) fprintf('%6.2f %8.4f %6.2f %13.4e %13.4e \n',tab') x_val=2; j=round((x_val-x0)/h)+1; tab1=[t sol1(:,j) sol2(:,j) err1(:,j) err2(:,j)]; tab1_rid=[tab1(1:2:end,:);tab1(end,:)]; fprintf(' \n\n Tabella per x=2 \n\n t \t\t sol1 \t\t sol2 \t\t err1 \t\t err2 \n') fprintf(' %7.3f %14.6e %14.6e %10.2e %10.2e \n', tab1_rid')

46. Istruzioni per la rappresentazione grafica h1=num2str(h); k1=num2str(k); titolo1=['metodo upwind-h=',h1,', k=', k1]; titolo2=['metodo implicito-h=',h1,', k=',k1]; figure(1) surf(x,t,Uvera),xlabel('x'),ylabel('t'),title('Soluzione vera') figure(2) surf(x,t,sol1),xlabel('x'),ylabel('t') title(['Soluzione approssimata:',titolo1]) figure(3) surf(x,t,sol2),xlabel('x'),ylabel('t') title(['Soluzione approssimata:',titolo2]) figure(4) contour(x,t,sol1,20),xlabel('x'),ylabel('t') title(['Curve di livello:',titolo1]) figure(5) contour(x,t,sol2,20),xlabel('x'),ylabel('t') title(['Curve di livello:',titolo2])

47. Tabelle dei risultati: quesiti 2) e 3) Tabella per x=2 t sol1 sol2 err1 err2 0.000 4.000000e+000 4.000000e+000 0.00e+000 0.00e+000 0.080 3.717258e+000 3.726409e+000 2.48e-002 3.39e-002 0.160 3.454842e+000 3.471142e+000 4.63e-002 6.26e-002 0.240 3.211372e+000 3.233120e+000 6.49e-002 8.66e-002 0.320 2.985532e+000 3.011265e+000 8.09e-002 1.07e-001 0.400 2.776066e+000 2.804522e+000 9.48e-002 1.23e-001 0.480 2.581778e+000 2.611861e+000 1.07e-001 1.37e-001 0.560 2.401536e+000 2.432298e+000 1.17e-001 1.47e-001 0.640 2.234266e+000 2.264889e+000 1.25e-001 1.56e-001 0.720 2.078954e+000 2.108743e+000 1.32e-001 1.62e-001 0.800 1.934646e+000 1.963017e+000 1.37e-001 1.66e-001 0.880 1.800444e+000 1.826920e+000 1.41e-001 1.68e-001 0.960 1.675506e+000 1.699703e+000 1.44e-001 1.68e-001 1.000 1.616261e+000 1.639202e+000 1.45e-001 1.68e-001 h k k+h err_imp err_upw 0.20 0.0400 0.24 3.6397e-001 2.4374e-001

48. Rappresentazione della soluzione vera

49. Superficie approssimata: metodo upwind

50. Curve di livello: metodo upwind