10a lezione di laboratorio

Laurea Ingegneria CIVILE

Lauree Specialistiche in Ingegneria CHIMICA, ELETTRONICA,

AMBIENTE

Laurea Ingegneria CIVILE

Lauree Specialistiche in Ingegneria CHIMICA, ELETTRONICA,

AMBIENTE

a.a. 2007-2008

Generalità su un problema del 2° ordine di tipo ellittico

),(),(),(

),(),(),(

yxyxgyxu

yxyxfyxu

),(),(),(

),(),(),(

yxyxgyxu

yxyxfyxu

.),(),,(:),( dcybaxyx

2

2

2

2

yx

Sia dato il seguente problema di Poisson- Dirichlet:

operatore laplaciano

con:

Approssimazioni utilizzatePartizioniamo con maglie quadrate di passo h, per cui a=Nh e b=Mh.

Nei nodi interni

si approssimano le derivate parziali seconde con le differenze finite:

,1,...,2,1,1,...,2,1per),( MjNiyx ji

l’errore locale è:

;2

,2

2

1,,1,

),(

2

2

2

,1,,1

),(

2

2

h

uuu

y

u

h

uuu

x

u jijiji

yx

jijiji

yx jiji

;

2,

22

1,,1,

),(

2

2

2

,1,,1

),(

2

2

h

uuu

y

u

h

uuu

x

u jijiji

yx

jijiji

yx jiji

).( 2hO

Costruzione del problema discreto2

, 1, 1, , 1 , 1 ,

0, 0, , ,

,0 ,0 , ,

4 1, 2, ..., 1 ; 1, 2, ..., 1

(*) , 0,1, ...,

, 0,1, ..., .

i j i j i j i j i j i j

j j N j N j

i i i M i M

u u u u u h f i N j M

u g u g j M

u g u g i N

dove fi,j = f( xi , yj ) e gi,j = g( xi , yj ), i=0,…,N, j=0,…,M.

Se si introduce il vettore:

si ottiene da (*), un sistema di

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.

1,1 2,1 1,1 1,2 2,2 1,2 1, 1 2, 1 1, 1, , ; , , ; ; , ,N N M M N Mu u u u u u u u u

)1()1( MN

Matrice A scritta a blocchi: caso N=5, M=4

,

0

0

CI

ICI

IC

A

1000

0100

0010

0001

,

4100

1410

0141

0014

IC

Risoluzione del sistemametodi 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

iterativi da utilizzarsi per matrici a blocchi.

Function PDE_ellittichefunction [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 commentog3=[g3,'+0*x']; g3=inline(g3);% vedi commentog4=[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 AR=[];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 soluziones=(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),-(N-% 1)); 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:

2 2

2 2, , 3, 3 3,3

( , ) 0, se 3 oppure se 3

u up x y x y

x y

u x y x y

2 2

2 2, , 3, 3 3,3

( , ) 0, se 3 oppure se 3

u up x y x y

x y

u x y x y

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 sol1 sol2

-3.000 0.00000000e+000 0.00000000e+000

-2.500 -2.44487048e-002 -2.54811702e-002

-2.000 -6.68102136e-002 -7.11188490e-002

-1.500 -2.84782305e-001 -2.93769126e-001

-1.000 -8.93486193e-001 -8.98795939e-001

-0.500 -1.63645642e+000 -1.63158849e+000

0.000 -1.98583980e+000 -1.98421091e+000

0.500 -2.13819801e+000 -2.13986488e+000

1.000 -2.31812451e+000 -2.31114343e+000

1.500 -2.13063817e+000 -2.12232331e+000

2.000 -1.50261774e+000 -1.49974215e+000

2.500 -7.41169824e-001 -7.40994870e-001

3.000 0.00000000e+000 0.00000000e+000

x sol1 sol2

-3.000 0.00000000e+000 0.00000000e+000

-2.500 -2.44487048e-002 -2.54811702e-002

-2.000 -6.68102136e-002 -7.11188490e-002

-1.500 -2.84782305e-001 -2.93769126e-001

-1.000 -8.93486193e-001 -8.98795939e-001

-0.500 -1.63645642e+000 -1.63158849e+000

0.000 -1.98583980e+000 -1.98421091e+000

0.500 -2.13819801e+000 -2.13986488e+000

1.000 -2.31812451e+000 -2.31114343e+000

1.500 -2.13063817e+000 -2.12232331e+000

2.000 -1.50261774e+000 -1.49974215e+000

2.500 -7.41169824e-001 -7.40994870e-001

3.000 0.00000000e+000 0.00000000e+000

figure(1):curve di livellopeaks(x,y), soluz._h=125

figure(2):soluz._h=0.25,soluz._h=0.125

Figure(3):peaks(x,y)

Funzione peaks(x,y)

%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

%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

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:

Rappresentazione della funzione peaks(x,y)

-3-2

-10

12

3

-2

0

2

-10

-5

0

5

x

3 (1-x).2 exp(-(x.2) - (y+1).2) -...- 1/3 exp(-(x+1).2 - y.2)

y

-6

-4

-2

0

2

4

6

8

Esercizio 2: parziale di Metodi Numerici per l’Ingegneria 30/11/2004

Si consideri il seguente problema alle derivate parziali:

1 -Si stabilisca, motivando la risposta, il tipo di equazione alle derivate parziali . Si verifichi inoltre, che la funzione

è soluzione del problema dato.

xyeyxU ),(xyeyxU ),(

]1,0[),1(,1),0(

]1,0[)1,(,1)0,(

)1,0(),1,0(0)( 22

yeyuyu

xexuxu

yxuyxu

y

x

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, 43 - 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: rappresentazionifigure(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: 5.97e-005

ERRORE MASSIMO PER h=1/8: 1.82e-005

___________________________________________________________

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: 5.97e-005

ERRORE MASSIMO PER h=1/8: 1.82e-005