10 a lezione di laboratorio Laurea Specialistica in Ingegneria Matematica Ingegneria dei Sistemi...

Post on 01-May-2015

214 views 0 download

Transcript of 10 a lezione di laboratorio Laurea Specialistica in Ingegneria Matematica Ingegneria dei Sistemi...

10a lezione di laboratorio

Laurea Specialistica in

Ingegneria MatematicaIngegneria dei Sistemi Energetici

Laurea Specialistica in

Ingegneria MatematicaIngegneria dei Sistemi Energetici

a.a. 2007-2008

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*picolorbar;title(' 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*picolorbar;title(' ezsurf')

-2-1

01

2

-2

-1

0

1

2-0.5

0

0.5

X

ezsurf

Y

-0.4

-0.3

-0.2

-0.1

0

0.1

0.2

0.3

0.4

Comando ezcontour

X

Y

X exp(-(X.2+Y.2))

-2 -1.5 -1 -0.5 0 0.5 1 1.5 2-2

-1.5

-1

-0.5

0

0.5

1

1.5

2

Con il comando contour tracciamo le linee di livello nel dominiofissato, 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]);

Comando ezsurf

Il comando ezsurf permette 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)

-2-1

01

2

-2-1

01

2

-5

0

5

x

x = 2 cos(s), y = 2 sin(s), z = z

y

z

% 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));

% 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);colorbartitle('surf')

% comando surf figure(4)surf(X,Y,Z);colorbartitle('surf')

-2-1

01

2

-2

-1

0

1

2-0.5

0

0.5

surf

-0.4

-0.3

-0.2

-0.1

0

0.1

0.2

0.3

0.4

Comando surf

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 curvetitle('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 curvetitle('contour')

contour

-2 -1.5 -1 -0.5 0 0.5 1 1.5 2-2

-1.5

-1

-0.5

0

0.5

1

1.5

2

Comando quiver

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

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

-2.5 -2 -1.5 -1 -0.5 0 0.5 1 1.5 2-2.5

-2

-1.5

-1

-0.5

0

0.5

1

1.5

2

2.5quiver

quiver e contour

-2 -1.5 -1 -0.5 0 0.5 1 1.5 2-2

-1.5

-1

-0.5

0

0.5

1

1.5

2

figure(5)contour(X,Y,Z,20;hold onquiver(X,Y,px,py);hold off

figure(5)contour(X,Y,Z,20;hold onquiver(X,Y,px,py);hold off

Esercizio 1

Sia dato il seguente problema alle derivate parziali (pde):

2

2

3 0 0,10 , 0

,0 2 exp 2 2

0, 3 2 exp 2 3 2

t xu u x t

u x x x x

u t t t - t

2

2

3 0 0,10 , 0

,0 2 exp 2 2

0, 3 2 exp 2 3 2

t xu u x t

u x x x x

u t t t - t

Quesiti a, b

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

a) Si verifichi che la funzione:

è soluzione del problema. 21 , 3 2 exp 2 3 2u x t x t x t 21 , 3 2 exp 2 3 2u x t x t x t

Soluzione del quesito a): Verifica

21 , 3 2 exp 2 3 2 ;u x t x t x t

2 2exp 2 3 2 1 4 3 2xu x t x t

2 23 exp 2 3 2 1 4 3 2tu x t x t

Quindi e:

3 0,t xu u

2,0 2 exp 2 2 ,u x x x

20, 3 2 exp 2 3 2 .u t t t

Quesito b): Metodo UPWINDApprossimazioni utilizzate:

, 1 , 1,

2

,0

2

0,

1 1,2, , , 0,1,

2 *exp 2* 2 0,1,2, ,

3 2 *exp 2* 3 2 0,1,

i j i j i j

i i i

j j j

u u u i N j

u x x i N

u t t j

30

ck k

h h

, 1 , , 1,

, ,,

i j i j

i j i j i j i j

x t x t

u u u uu u

t k x h

Indicando quindi si ottiene:

Sistema relativo al problema discretoSe si assume per ogni livello temporale j:

1j j jU AU V

1, 2, ,, , , ,Tj j j N jU u u u

0, ,0, ,0 ,0, ,0 ,T Tj j jV u t

Sappiamo che se:

1,ck

h

il problema discreto diventa:

il metodo CONVERGE quando ., 0h k

Costruzione delle formuleDalle relazioni:

1 1 1

2 1 2

1

1 0

0.

0 1 0

j j j

j j

N j N j

u u

u u

u u

si ottiene, per ogni livello temporale j, tenendo anche conto della condizione al contorno, il sistema lineare:

, 1 , 1,1 1,2, , , 0,1,i j i j i ju u u i N j

Forma della matrice A

1 0

,

0 1

A

Si deduce allora che la matrice A, di dimensioni NxN ,per ogni j, ha la forma:

e, nell’ipotesi risulta:

11 1 1.A A

0 1,

Il metodo è condizionatamente stabile!!!

Quesito b): Metodo IMPLICITO

In questo caso si colloca la pde in

3,

ck k

h h

il problema diventa:

1 1

, 1 , , 1 1, 1

, ,

,i j i j

i j i j i j i j

x t x t

u u u uu u

t k x h

Si indica ancora:

1,i jx t

e si approssimano le derivate parziali con:

due differenze all’indietro!!

Problema discretizzato per il metodo implicito

, 1 1, 1 ,

2

,0

2

0,

1 1,2, , , 0,1,

2 *exp 2* 2 0,1, ,

3 2 *exp 2 3 2 0,1,

i j i j i j

i i i

j j j

u u u i N j

u x x i N

u t t j

Indicando ancora: 1, 2, ,, , , ,Tj j j N jU u u u

1 0, 1 1,0, ,0 ,0, ,0 ,T

Tj j jV u t

si ottiene il sistema:

Sistema relativo al metodo implicito

1 1 sistema lineare,j j jAU U V 1 0

,

0 1

A

Poiché risulta:

11 j jE A E

La formula per l’errore è:

, !!!h kil metodo converge quindi per:

1

11 , ,A h k

Convergenza incondizionata!!!

Istruzioni relative al quesito b)clear all; clct0=0;x0=0;xN=10;h=0.25; M=10;c='3';c1=eval(c);f='(x-2).*exp(-2*(x-2).^2)';% cond. inizialeg='-(3*t+2).*exp(-2*(3*t+2).^2)';%cond.contornor='0';%termine notofprintf(['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

Function PDE_upwindx=(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;endt=linspace(t0,tM,M+1)';x=[x0;x];

Function PDE_implicitox=(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;endt=linspace(t0,tM,M+1)';x=[x0;x];

Risultati al variare del passo k

h k Errore max IMP Err max UPW h+k 0.25 0.05 0.6 0.2467 0.1249 0.30 0.25 h/3 1 0.2878 4.4409e-016 0.3333 0.25 0.1 1.2 0.2957 0.3387 0.35 0.25 0.5 6 0.3523 9.9564e+006 0.75

l’ implicito converge, upwind è instabile! e quindi non converge.

1:

entrambi i metodi sono consistenti:

sono stabili, e quindi convergono.

1: O h k

M = 10

Osservazione sul caso

1 11 = coeff. angolare delle caratteristiche!!!

3

k

h c

1

il metodo upwind fornisce:

1j j jU AU V

0 0

1con .

0 1 0

A

03x t x Linee caratteristiche :

Commenti sul caso

1 1, 2, 1,

1 1,0 2,0 1,0

se 0 0

Tj j j j N j

T

N

U t u u u

j U u u u

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

1 Si ottiene, per la forma di ,A

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); %surfltitle('soluzione approssimata:metodo upwind')xlabel('x'),ylabel('t')

figure(2)C=contour(X,T,sol1,20); %20 curve di livellotitle('Curve di livello') xlabel('x'),ylabel('t')

Superficie: metodo upwind

0

5

10

0

0.5

1-0.4

-0.2

0

0.2

0.4

x

Soluzione approssimata: metodo upwind

t

h=0.25, k=h/3h=0.25, k=h/3

2 4 6 8 10-0.4

-0.3

-0.2

-0.1

0

0.1

0.2

0.3

0.4

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.

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

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

Esercizio 2

Si consideri il seguente problema misto ai valori

iniziali ed al contorno, con coefficienti non costanti:

2

4

2 0 0,3 , 0

1( ,0) 0

11

(0, ) 01

t xu tu x t

u x f x xx

u t g t tt

Quesiti a, b

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.

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

Soluzione del quesito a)

2 20 0 ,x t x x x t p x t

0, con ,0x t x

02 0dx

t x xdt

Si verifichi ora che

è la soluzione del problema pde+ condizione iniziale che soddisfa anche la condizione assegnata al bordo.

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

222

1,

1u x t f x t

x t

Soluzione del quesito b)clear all; clct0=0;tM=3; % in questo caso si assegna tmaxx0=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 veratab=[];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];endtab=[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')

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

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 livelloset(gca, 'FontWeight','bold','Fontsize',12)set(H,'LineWidth',2)title(['Curve di livello:',titolo])xlabel('x'); ylabel('t')

Rappresentazione della soluzione vera

Superficie approssimata: metodo upwind

k=h/6k=h/6

Curve di livello: metodo upwind

k=h/6k=h/6

Superficie approssimata: metodo implicito

k=h/6k=h/6

Curve di livello: metodo implicito

k=h/6k=h/6

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

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 quindi k=h/6.

Esercizio 3

Sia dato il seguente problema alle derivate

parziali a coefficienti non costanti:

.00),0(

,40)0,(

0,40)2()(2

22

ttu

xxxu

txxtxeuxtu txt

.00),0(

,40)0,(

0,40)2()(2

22

ttu

xxxu

txxtxeuxtu txt

con soluzione vera:

.),()1( 2xetxu t .),()1( 2xetxu t

Quesiti 1) e 2)

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

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

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.

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

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

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

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

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

h k k+h err_imp err_upw

0.20 0.0400 0.24 3.6397e-001 2.4374e-001

Rappresentazione della soluzione vera

Superficie approssimata: metodo upwind

Curve di livello: metodo upwind

Superficie approssimata: metodo implicito

Curve di livello: metodo implicito