aa
7 8
lezione di laboratorio
Laurea Ingegneria CIVILE
Lauree Specialistiche in Ingegneria
CHIMICA, ELETTRONICA,
AMBIENTE
a.a. 2007-2008
Comando ezsurf
ezsurf
0.4
0.3
0.5
0.2
0.1
0
0
-0.1
-0.5
2
-0.2
1
2
0
1
-0.3
0
figure(1)
-0.4
-1
-1
z='X.*exp(-(X.^2+Y.^2))'; Y
-2 -2
X
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')
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))'; 2
ezcontour(z,[-2,2,-2,2]);
X exp(-(X. 2+Y.2))
1.5
1
Y
0.5
0
-0.5
-1
-1.5
-2
-2
-1.5
-1
-0.5
0
0.5
1
1.5
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)
x = 2 cos(s), y = 2 sin(s), z = z
z
5
0
-5
2
1
1
0
-1
y
-2
-2
0
-1
x
2
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));
surf
0.4
0.3
0.5
% comando surf
figure(4)
surf(X,Y,Z);colorbar
0
title('surf')
0.2
0.1
0
-0.1
-0.5
2
-0.2
1
2
1
0
0
-1
-1
-2
-2
-0.3
-0.4
Comando contour
contour
2
1.5
1
0.5
0
-0.5
-1
-1.5
-2
-2
-1.5
-1
figure(5)
contour(X,Y,Z,20)
-0.5
0
0.5
1
1.5
2
% 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')
Comando quiver
quiver
2.5
2
1.5
1
0.5
0
-0.5
-1
-1.5
-2
-2.5
-2.5
-2
-1.5
-1
-0.5
0
0.5
1
1.5
2
figure(4)
[px,py]=gradient(Z,.2,.2);%[px,py]=gradient(Z);
quiver(X,Y,px,py)
title('quiver')
quiver e contour
2
1.5
1
0.5
0
-0.5
-1
-1.5
-2
-2
-1.5
-1
-0.5
0
figure(5)
contour(X,Y,Z,20;hold on
quiver(X,Y,px,py);hold off
0.5
1
1.5
2
Esercizio 1
Sia dato il seguente problema alle derivate
parziali (pde):
 ut  3ux  0 x   0,10  , t  0

 u  x , 0     x    x  2  exp 2  x  2  2


2
 u  0, t     t     3t  2  exp -2  3t  2 





Quesiti a, b
a) Si verifichi che la funzione:
2
1 u  x, t    x  3t  2 exp 2  x  3t  2


è soluzione del problema.
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.
Soluzione del quesito a): Verifica
1 u  x, t    x  3t  2 exp  2  x  3t  2

ux  exp 2  x  3t  2

2
;
 1  4  x  3t  2 
ut   3 exp 2  x  3t  2 
Quindi
2
2
2
 1  4  x  3t  2 
2
ut  3ux  0, e:


u  x,0   x  2 exp 2  x  2 ,
2


u  0, t     3t  2 exp 2  3t  2 .
2
Quesito b): Metodo UPWIND
Approssimazioni utilizzate:
ui , j 1  ui , j
u

,
x
,
t
t  i j 
k
ui , j  ui 1, j
u

x
,
t
x  i j 
h
ck 3k

 0 si ottiene:
Indicando quindi  
h
h

u i , j 1  1    u i , j   u i 1, j


2
 u i ,0   xi  2 *exp 2*  xi  2 

 u   3t  2 *exp 2* 3t  2
j
j
 0, j





i  1, 2, , N , j  0,1,

i  0,1, 2, , N

2
j  0,1,
Problema discreto
Se si assume per ogni livello temporale j:
T


U j   u1, j , u2, j , , uN , j  ,
V j   u0, j ,0,
 
T

,0   t j ,0,

T

,0 ,

il problema discreto diventa:
U j 1  AU j   V j
Sappiamo che se:   ck  1,
h
il metodo CONVERGE quando h, k  0 .
Costruzione delle formule
Dalle relazioni:
ui , j 1  1    ui , j   ui 1, j
i  1,2, , N , j  0,1,
si ottiene, per ogni livello temporale j,
tenendo anche conto della condizione al
contorno la seguente equazione vettoriale:
 u1 j 1  1  

 
 u2 j 1   



 
 uN j 1   0

 j 
0   u1 j 
 
u 
  2j    0 .

 


 

1     uN j 
 0 
Forma della matrice A
Si deduce allora che la matrice A, di
dimensioni NxN ,per ogni j, ha la forma:
1  
 
A


 0
e, nell’ipotesi

0 

,


1 
0    1, risulta:
A   A 1  1      1      1.
Il metodo è condizionatamente stabile!!!
Quesito b): Metodo IMPLICITO
In questo caso si colloca la pde in
 x ,t 
i
j 1
e si approssimano le derivate parziali con:
ui , j 1  ui , j
u

,
t  xi ,t j1 
k
ui , j 1  ui 1, j 1
u

x  xi ,t j1 
h
due differenze all’indietro!!
ck 3k
Si indica ancora:   h  h ,
il problema diventa:
Problema discretizzato con il metodo
implicito

1    ui , j 1   ui 1, j 1  ui , j


2
 ui ,0   xi  2  * exp 2*  xi  2 

 u   3t  2 * exp 2 3t  2 2
j
j
 0, j




i  1, 2,


i  0,1,

Indicando ancora: U j   u1, j , u2, j ,
V j 1   u0, j 1 , 0,
 
, 0   t j 1 , 0,
T
, N , j  0,1,
,N
j  0,1,
T

, uN , j  ,
T
, 0 ,
si ottiene il sistema:
Sistema relativo al metodo implicito
AU j 1  U j  V j 1  sistema lineare,
1  
 
A


 0

Poiché risulta:
0 

,
 La formula per l’errore è:

1
1 
E j 1  A E j
1
A
1
 1 h, k ,
il metodo converge quindi per:
h, k !!!
Convergenza incondizionata!!!
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
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];
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];
Risultati al variare del passo k
h
0.25
0.25
0.25
0.25
k
0.05
h/3
0.1
0.5

0.6
1
1.2
6
M = 10
Errore max IMP Err max UPW
0.2467
0.1249
0.2878
4.4409e-016
0.2957
0.3387
0.3523
9.9564e+006
h+k
0.30
0.3333
0.35
0.75
  1: entrambi i metodi sono consistenti:  O  h  k 
sono stabili, e quindi convergono.
  1: l’ implicito converge, upwind è instabile!
e quindi non converge.
 1
Osservazione sul caso
Linee caratteristiche : x  3t  x0
k 1 1
  1    = coeff. angolare delle caratteristiche!!!
h 3 c
il metodo upwind fornisce:
0
1
con A  


0
1
0

.


0
U j 1  AU j  Vj
Commenti sul caso   1
Si ottiene, per la forma di
A,
U j 1
 
  t j
u1, j
T

uN 1, j 
u2, j

se j  0 U1    0  u1,0
u2,0
uN 1,0 
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!!!!
T
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')
Superficie: metodo upwind
h=0.25, k=h/3
Soluzione approssimata: metodo upwind
0.4
0.2
0
-0.2
-0.4
1
10
0.5
t
5
0 0
x
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.
0.4
0.3
0.2
0.1
0
-0.1
-0.2
-0.3
-0.4
2
4
6
8
10
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
Esercizio 2
Si consideri il seguente problema misto ai valori
iniziali ed al contorno, con coefficienti non costanti:

 ut  2tux  0

1

 u( x , 0)  f  x  
2
1 x

1

 u(0, t )  g  t   1  t 4
x   0, 3 , t  0
x0
t0
Quesiti a, b
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.
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.
Soluzione del quesito a)
Per individuare le caratteristiche della pde data, si risolve
il problema di Cauchy:
 dx
  2t x  0   x0
 dt
Esso ha la soluzione
x  t 2  x0  x0  x  t 2  p  x, t 
 x, t 
e p(x,t) è la linea che collega
Si verifichi ora che

u  x, t   f x  t
2
con  x0 ,0

1


1 x  t
è la soluzione del problema pde+ condizione iniziale che
soddisfa anche la condizione assegnata al bordo.
2
2
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')
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 livello
set(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/6
Curve di livello: metodo upwind
k=h/6
Superficie approssimata: metodo
implicito
k=h/6
Curve di livello: metodo implicito
k=h/6
Errori in t= tM=3
h
k
k+h alpha err_imp
0.20 0.0333 0.2333 1.00 3.1685e-001
0.10 0.0167 0.1167 1.00 2.1497e-001
0.05 0.0083 0.0583 1.00 1.3615e-001
err_upw
2.0392e-001
1.2541e-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:
 ut  ( t 2  x )u x  e  t ( x  2t 2 ) x 0  x  4, t  0

2
u
(
x
,
0
)

x
0  x  4,

 u(0, t )  0
t  0.

con soluzione vera:
t
(1) u( x , t )  e x .
2
Quesiti 1) e 2)
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.
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”.
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
0.000
4.000000e+000
0.080
3.717258e+000
0.160
3.454842e+000
0.240
3.211372e+000
0.320
2.985532e+000
0.400
2.776066e+000
0.480
2.581778e+000
0.560
2.401536e+000
0.640
2.234266e+000
0.720
2.078954e+000
0.800
1.934646e+000
0.880
1.800444e+000
0.960
1.675506e+000
1.000
1.616261e+000
h
sol2
4.000000e+000
3.726409e+000
3.471142e+000
3.233120e+000
3.011265e+000
2.804522e+000
2.611861e+000
2.432298e+000
2.264889e+000
2.108743e+000
1.963017e+000
1.826920e+000
1.699703e+000
1.639202e+000
k
0.20 0.0400
k+h
0.24
err1
0.00e+000
2.48e-002
4.63e-002
6.49e-002
8.09e-002
9.48e-002
1.07e-001
1.17e-001
1.25e-001
1.32e-001
1.37e-001
1.41e-001
1.44e-001
1.45e-001
err_imp
3.6397e-001
err2
0.00e+000
3.39e-002
6.26e-002
8.66e-002
1.07e-001
1.23e-001
1.37e-001
1.47e-001
1.56e-001
1.62e-001
1.66e-001
1.68e-001
1.68e-001
1.68e-001
err_upw
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