a
10
lezione di laboratorio
Laurea Ingegneria CIVILE
Lauree Specialistiche in Ingegneria
CHIMICA, ELETTRONICA,
AMBIENTE
a.a. 2007-2008
Generalità su un problema del 2°
ordine di tipo ellittico
Sia dato il seguente problema di Poisson- Dirichlet:
  u( x , y )  f ( x , y ) ( x , y )  

 u( x , y )  g ( x , y ) ( x , y )  
con:


 2 2
x
y
2
2
operatore laplaciano
  ( x, y) : x  (a, b), y  (c, d ).
Approssimazioni utilizzate
Partizioniamo    con maglie quadrate di
passo h, per cui a=Nh e b=Mh.
Nei nodi interni
( x i , y j ) per i  1,2,..., N  1,
j  1,2,..., M  1,
si approssimano le derivate parziali seconde
con le differenze finite:
ui 1, j  2ui , j  ui 1, j
ui , j 1  2ui , j  ui , j 1
 2u
 2u

,

;
2
2
2
2
x ( x , y )
h
y ( x , y )
h
i
j
i
2
l’errore locale è: O( h ).
j
Costruzione del problema discreto
4 ui , j  ui 1, j  ui 1, j  ui , j 1  ui , j 1  h2 f i , j i  1, 2,..., N  1; j  1, 2,..., M  1

(*) u0, j  g0, j , uN , j  gN , j
j  0,1,..., M

i  0,1,..., N .
ui ,0  gi ,0 , ui , M  gi , M
dove fi,j = f( xi , yj ) e gi,j = g( xi , yj ), i=0,…,N, j=0,…,M.
Se si introduce il vettore:
u1,1 , u2,1 ,
uN 1,1 ; u1,2 , u2,2 ,
uN 1,2 ;
; u1, M 1 , u2, M 1 ,
uN 1, M 1
si ottiene da (*), un sistema di ( N  1)  ( M  1)
equazioni in altrettante incognite con matrice A dei
coefficienti simmetrica e a dominanza diagonale debole.
Si può verificare che A è definita positiva, dunque il
sistema e’ univocamente risolubile.
Matrice A scritta a blocchi:
caso N=5, M=4
C

A   I
 0
I
C
I
0 

 I ,
C 
0
 4 1 0
1 0
 1 4  1 0 
0 1
, I  
C
 0  1 4  1
0 0



0 1 4 
0
0 0
0 0
0 0
1 0

0 1
Risoluzione del sistema
metodi idonei per la risoluzione del sistema
tenendo conto delle caratteristiche della matrice A:
a) se A è di piccole dimensioni:
1) Cholesky,
2) Gauss senza pivoting;
b) se A è di grandi dimensioni:
1) Gauss-Seidel,
2) SOR (da preferire perché si conosce il parametro
ottimale).
Sono state costruite versioni degli algoritmi dei metodi
Function PDE_ellittiche
function [x,y,sol]=PDE_ellittiche(y0,yM,x0,xN,h,f,r,
g1,g2,g3,g4);
x=(x0:h:xN)'; x(end)=xN; N=length(x)-1;
y=(y0:h:yM)'; y(end)=yM; M=length(y)-1;
f=[f,'+0*x+0*y']; f=inline(f);% vedi commento
r=[r,'+0*x+0*y']; r=inline(r);% vedi commento
g1=[g1,'+0*x']; g1=inline(g1);% vedi commento
g2=[g2,'+0*y']; g2=inline(g2);% vedi commento
g3=[g3,'+0*x']; g3=inline(g3);% vedi commento
g4=[g4,'+0*y']; g4=inline(g4);% vedi commento
% Queste istruzioni sono utilizzate solo se le funzioni f,
% r,g1,g2,g3,g4 sono costanti.
% Costruzione della matrice dei coefficienti A
R=[];
for j=1:M-1
R=[R,feval(r,x(2:N),y(j+1))'];
end
% s è il numero di nodi interni in cui viene calcolata la
soluzione
s=(M-1)*(N-1);
Function PDE_ellittiche
%
%
%
%
%
%
%
N.B. Le righe commentate possono utilizzarsi in caso di matrici
di piccole dimensioni. In caso di grandi dimensioni conveniente
usare la costruzione mediante il comando spdiags.
A=diag((4-h^2*R))-diag(ones(s-1,1),1)-diag(ones(s-1,1),-1);
A=A-diag(ones((M-2)*(N-1),1),N-1)-diag(ones((M-2)*(N-1),1),-(N1));
e = ones(s,1);
A = spdiags([ -e -e (4-h^2*R)' -e -e],[-(N-1) -1 0 1 N-1],s,s);
% devono essere azzerati alcuni elementi delle codiagonali
% principali
for k=N-1:N-1:(M-2)*(N-1)
A(k,k+1)=0;
A(k+1,k)=0;
end
% costruzione del vettore termini noti
F=[];
for j=1:M-1
F=[F,feval(f,x(2:N),y(j+1))'];
end
Function PDE_ellittiche
% condizioni al bordo inferiore
G=zeros(1,(N-1)*(M-1));
G(1:N-1)=[feval(g1,x(2:N))'];
G(1)=G(1)+feval(g4,y(2));
G(N-1)=G(N-1)+feval(g2,y(2));
% condizioni ai bordi destro e sinistro
for j=2:M-2
G((j-1)*(N-1)+1)=feval(g4,y(j+1));
G(j*(N-1))=feval(g2,y(j+1));
end
% condizioni al bordo superiore
G((M-2)*(N-1)+1:end)=[feval(g3,x(2:N))'];
G((M-2)*(N-1)+1)=G((M-2)*(N-1)+1)+feval(g4,y(M));
G((M-1)*(N-1))=G((M-1)*(N-1))+feval(g2,y(M));
F=h^2*F+G;
Function PDE_ellittiche
% calcolo della soluzione
U=A\F';
%la soluzione trovata viene messa in una matrice
for j=1:M-1
sol(j,:)=[feval(g4,y(j+1))
U((j-1)*(N-1)+1:(j)*(N-1))' feval(g2,y(j+1))];
end
sol=[feval(g1,x)';sol;feval(g3,x)' ];
Esercizio 1:esame di Metodi
Numerici per l’Ingegneria (19/12/2005)
Sia p(x,y) la funzione definita dal file di Matlab
peaks(x,y)e si consideri il seguente problema alle
derivate parziali:
 u  u
 2  2  p  x , y   x , y    3, 3   3, 3
y
 x
 u( x , y )  0, se x  3 oppure se y  3

2
2
1) Determinare, motivando la risposta, a quale classe
appartiene il problema proposto.
Quesito 2)
2) Si costruisca un file MATLAB:
Cognome_nome_studente.m che, una volta avviato:
a) faccia visualizzare una schermata con i dati personali
e una breve presentazione del problema;
b) permetta di dare in input i passi lungo l’asse x e y
dati da: h1=1/4, h2=1/8 ;
c) determini la soluzione utilizzando il metodo alle
differenze finite;
d) costruisca una tabella relativa al primo ed al secondo
caso, che riporti intestazione:
x sol1 sol2
Quesito 2)
nella prima colonna i valori , comuni nei due casi ,
nelle successive colonne i valori corrispondenti della
soluzione approssimata per y = 0; si stampi la tabella ogni
2 righe e si utilizzino i seguenti formati di stampa:
3 cifre decimali e formato virgola fissa per i valori dei
nodi ,
8 cifre decimali e formato floating point per le soluzioni
approssimate nei due casi.
Quesiti 3) e 4)
3) Si eseguano due figure (relative ai due casi) in cui si
utilizzi il comando subplot con due finestre grafiche:
nelle due finestre della prima figura, si riportino
rispettivamente, 20 linee di livello di p(x,y) e di sol2;
nella seconda figura si riportino le rappresentazioni di
sol1 e sol2 rispettivamente nella prima e seconda finestra.
In una terza figura si riporti la rappresentazione della
soluzione vera.
4) Si commentino i risultati specificando, in particolare,
se essi soddisfano le aspettative teoriche.
Istruzioni Quesito 2)
clc; clear all
x0=-3;xN=3; y0=-3;yM=3;
h=[1/4 1/8];
f='-peaks(x,y)';%attenzione al segno
r='0';g1='0';g2='0';g3='0';g4='0'; y_val=0;
tab=[]; figure(2)
for i=1:length(h)
[x,y,sol]=PDE_ellittiche(y0,yM,x0,xN,h(i),f,r,g1,g2,g3,g4);
step=round(h(1)/h(i));
ind_val=round((y_val-y0)/h(i))+1;
tab=[tab sol(ind_val,1:step:end)'];
subplot(2,1,i),surf(x,y,sol);title(['sol.–h= ' num2str(h(i))]);
xlabel('x');ylabel('y');axis([-3 3 -3 3 -5 1])
end
tab=[x(1:step:end) tab];
fprintf('
x
sol1
sol2 \n')
fprintf(' %8.3f
%16.8e
%16.8e \n',tab(1:2:end,:)')
Istruzioni Quesito 3): curve di livello
[x,y]=meshgrid(x,y);
p=-eval(f);
figure(1)
subplot(2,1,1),contour(x,y,p,20)
title('peaks(x,y)')
xlabel('x');ylabel('y')
axis([-3 3 -3 3 ])
subplot(2,1,2),contour(X,Y,sol,20)
title('sol.- h =',num2str(h(end)))
xlabel('x');ylabel('y')
axis([-3 3 -3 3 ])
Risultati: tabella per y=0
x
-3.000
-2.500
-2.000
-1.500
-1.000
-0.500
0.000
0.500
1.000
1.500
2.000
2.500
3.000
sol1
0.00000000e+000
-2.44487048e-002
-6.68102136e-002
-2.84782305e-001
-8.93486193e-001
-1.63645642e+000
-1.98583980e+000
-2.13819801e+000
-2.31812451e+000
-2.13063817e+000
-1.50261774e+000
-7.41169824e-001
0.00000000e+000
sol2
0.00000000e+000
-2.54811702e-002
-7.11188490e-002
-2.93769126e-001
-8.98795939e-001
-1.63158849e+000
-1.98421091e+000
-2.13986488e+000
-2.31114343e+000
-2.12232331e+000
-1.49974215e+000
-7.40994870e-001
0.00000000e+000
figure(1):curve di livello
peaks(x,y), soluz._h=125
figure(2):
soluz._h=0.25,soluz._h=0.125
Figure(3):peaks(x,y)
Funzione peaks(x,y)
Digitando help peaks in Matlab, appare una
descrizione di tale funzione. Digitando edit peaks,
appare il file che definisce peaks(x,y).
La rappresentazione di peaks(x,y) può quindi essere
ottenuta utilizzando le seguenti istruzioni:
%Definizione e rappresentazione di peaks(x,y)
>> z = '3*(1-x).^2.*exp(-(x.^2) - (y+1).^2)...
- 10*(x/5 - x.^3 - y.^5).*exp(-x.^2-y.^2)...
- 1/3*exp(-(x+1).^2 - y.^2)';
>> ezsurf(z,[-3 3 -3 3])
>> colorbar
Rappresentazione della funzione
peaks(x,y)
3 (1-x).2 exp(-(x.2) - (y+1).2) -...- 1/3 exp(-(x+1).2 - y.2)
8
6
5
4
0
2
0
-5
-2
-10
2
0
-2
y
-3
-2
0
-1
x
1
2
3
-4
-6
Esercizio 2: parziale di Metodi Numerici
per l’Ingegneria 30/11/2004
Si consideri il seguente problema alle derivate parziali:
  u  ( x 2  y 2 )u  0

x
 u( x ,0)  1, u( x , 1)  e
 u(0, y )  1, u(1, y )  e  y

x  (0,1),
x  [0,1]
y  [0,1]
y  (0,1)
1 -Si stabilisca, motivando la risposta, il tipo di equazione
alle derivate parziali . Si verifichi inoltre, che la
funzione
U ( x, y)  e
è soluzione del problema dato.
 xy
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;
• utilizzando il metodo alle differenze finite, con stessa spaziatura
lungo x e lungo y, calcoli, sui punti dei due reticoli corrispondenti
ai passi h1=1/4 e h2=1/8, la soluzione approssimata e l’errore
assoluto;
• faccia visualizzare, nei due casi, la matrice della soluzione
approssimata nei punti coincidenti per i due passi e l’ errore
assoluto massimo, con i seguenti formati di stampa:
6 cifre decimali e formato virgola fissa per la soluzione;
2 cifre decimali e formato esponenziale per l’errore assoluto
massimo.
Quesiti 3, 4
3 - In una figura, si riporti il grafico della soluzione
vera. In una seconda figura con 4 finestre grafiche si
riportino, nella prima e nella seconda finestra, le
soluzioni approssimate nei due casi e, nella terza e
nella quarta, gli errori con i due passi richiesti. In una
terza figura, con uno schema grafico uguale a quello
della figura precedente, si grafichino 20 curve di
livello relative alla soluzione approssimata nei casi
h1=1/4 e h2=1/8 e si riportino le curve di livello
dell’errore nei due casi. Si corredi ogni singola figura
con intestazione contenente la specifica del passo.
4 - Si commentino i risultati verificando se essi
soddisfano le aspettative teoriche.
Istruzioni quesito 2
clc; clear all
x0=0; xN=1; y0=0; yM=1;
h=[1/4 1/8];
f='0'; r='-x.^2-y.^2';
g1='1'; g2='exp(-y)'; g3='exp(-x)'; g4='1';
s='___________________________________________________________';
for i=1:length(h)
step=round(h(1)/h(i));
[x,y,sol]=PDE_ellittiche(y0,yM,x0,xN,h(i),f,r,g1,g2,g3,g4);
[X,Y]=meshgrid(x,y);
Uvera=exp(-X.*Y);
err=abs(Uvera-sol);errmax=max(max(err));
xr=x(1:step:end);yr=y(1:step:end);
solr=sol(1:step:end,1:step:end);
end
Istruzioni quesito 2: tabelle
(all’interno del ciclo for)
n=length(xr);str=[];
for l=1:n
str=[str,'%10.6f '];
end
str=[str,'\n'];
disp(s)
fprintf(['SOL. APPROSSIMATA CON h=' num2str(h(i)) ' \n '])
disp(s)
fprintf(str,solr')
fprintf('\n\n')
disp(s)
fprintf(['ERRORE MASSIMO PER h=' num2str(h(i)) ': %10.2e \n '],
errmax)
Istruzioni quesito 3: rappresentazioni
figure(2) %ISTRUZIONI ALL’INTERNO DEL CICLO FOR
subplot(2,2,i),surf(x,y,sol)
set(gca,'FontWeight','bold','Fontsize',12)
title(['Sol. approssimata - h= ' num2str(h(i)) ])
subplot(2,2,i+2), surf(x,y,err)
set(gca,'FontWeight','bold','Fontsize',12)
title(['Err - h= ' num2str(h(i)) ])
figure(3)
subplot(2,2,i), contour(x,y,sol,20)
set(gca,'FontWeight','bold','Fontsize',12)
title(['C. livello (sol)- h=' num2str(h(i)) ]),colorbar
subplot(2,2,i+2), contour(x,y,err,20)
set(gca,'FontWeight','bold','Fontsize',12)
title(['C. livello (err) - h=' num2str(h(i)) ]),colorbar
figure(1)
surf(X,Y,Uvera)
set(gca,'FontWeight','bold','Fontsize',12)
title('Soluzione vera: U(x,y)=exp(-xy)')
Risultati quesito 2
___________________________________________________________
SOLUZIONE APPROSSIMATA CON h=1/4
___________________________________________________________
1.000000
1.000000
1.000000
1.000000
1.000000
1.000000
0.939430
0.882531
0.829074
0.778801
1.000000
0.882531
0.778855
0.687349
0.606531
1.000000
0.829074
0.687349
0.569841
0.472367
1.000000
0.778801
0.606531
0.472367
0.367879
___________________________________________________________
SOLUZIONE APPROSSIMATA CON h=1/8
___________________________________________________________
1.000000
1.000000
1.000000
1.000000
1.000000
1.000000
0.939418
0.882507
0.829042
0.778801
1.000000
0.882507
0.778816
0.687306
0.606531
1.000000
0.829042
0.687306
0.569799
0.472367
1.000000
0.778801
0.606531
0.472367
0.367879
ERRORE MASSIMO PER h=1/4:
ERRORE MASSIMO PER h=1/8:
5.97e-005
1.82e-005