Esercitazione esame metodi numerici

Formattazione del testo con tag HTML

%VALORI NODALIdir=load('mesh0.bound');ndir=length(dir);

%DATI FUNZIONE ED EVENTUALE SOLUZIONE ALLO STAZIONARIOfun=@(x,y) (4-2*(x.^2)-2*(y.^2));ustaz=@(x,y) -(x.*y).^2+x.^2+y.^2-1; %soluzione allo stazionarioureal=ustaz(x,y);

%DATI PER IL GCMkmax=300;ncr=0;iprec=2;deltat=0.1;teta=1;tmax=3;ndt=1;tempo=0;S=1;end

3) TOPOLOGIA MATRICE RIGIDEZZASuccessivamente viene calcolata la topologia delle matrici H e P le quali saranno salvate in forma sparsaperché per griglie molto raffinate il salvataggio in forma sparsa ottimizza l’utilizzo della memoria e lavelocità del codice.function [ H,P ] = topol( elem,ne )%function con la quale calcolo la topologia dellamatrice di rigidezza e%capacità del miom problemarow=reshape(repmat(1:ne,3,1),1,[])';col=reshape(elem',1,[])';A=sparse(row,col,0);H=A'*A;P=H;end

4) ASSEMBLAGGIO MATRICE DI RIGIDEZZA E CAPACITA’Salvata la topologia delle matrici, posso calcolarmi i valori e assemblarli assieme. Il

Il termine generico delle matrici locali di rigidezza e capacità è: 4dove ∆ è l'area del singolo elemento finito e b e c sono valori ricavati dal calcolo delle funzioni base triangolari, le quali valgono 1 nel punto di appoggio e 0 negli altri punti:

Una volta calcolate le componenti locali, posso inserirle nella matrice globale di rigidezza elemento per elemento.

Il vettore del termine noto viene calcolato genericamente dalla formula:

L'integrale può essere anche ristretto all'unione dei triangoli che hanno in comune il vertice i; essendo baricentrico, possiamo assumere f(x,y) costante e uguale a f(xi,yi). La formula diventa:

La function implementata risulta essere:

function [ D,H,P,areanod,f ] = cicloelem( ne,nn,elem,x,y,S,H,P,fun )
%function per calcolare la matrice di rigidezza H e per P
%si usa un'implementazione agli elementi finiti triangolari
D=zeros(ne,1);
areanod=zeros(nn,1);
for k=1:ne
    e=elem(k,:);
    ii=e(1,1);
    jj=e(1,2);
    mm=e(1,3);
    % Calcolo di

b,b=zeros(3,1);c=zeros(3,1);b(1)=y(jj)-y(mm);b(2)=y(mm)-y(ii);b(3)=y(ii)-y(jj);c(1)=x(mm)-x(jj);c(2)=x(ii)-x(mm);c(3)=x(jj)-x(ii);
C=[1 x(ii) y(ii); 1 x(jj) y(jj); 1 x(mm) y(mm)];D(k)=det(C)/2;
%assemblaggio matrici locali
%input funzione,areanod
%Calcolo area elemento
%Tipo di nodi
for ii=1:3
i=elem(k,ii);
for jj=1:3
j=elem(k,jj);
Hloc=(b(ii)*b(jj)+c(ii)*c(jj))/(4*D(k));
Ploc=(D(k)/12)*[2 1 1; 1 2 1; 1 1 2]*S;
H(i,j)=H(i,j)+Hloc;
P(i,j)=P(i,j)+Ploc(ii,jj);
end
areanod(i)=areanod(i)+D(k)/3;
end
end f=-(fun(x,y).*areanod);
end 55) SOLUZIONE ALLO STAZIONARIO
Ora si può iniziare a risolvere il problema. In questo caso possiamo verificare che la soluzione allostazionario sia corretta implementando il metodo del gradiente coniugato modificato(GCM) e confrontandola soluzione numerica nei vari punti con la soluzione reale analitica.
function [ x,k,T,Tvero,rvero,r2 ] = gradconmod( A,b,E,kmax,ncr,m,iprec )
%function del gradiente coniugato modificato con implementati i
%precondizionatori di Jacobi e delladecomposta incompleta di Cholesky; con possibilità di inizializzare la soluzione con le correzioni residue. ```html

r=b;x=zeros(m,1);T=zeros(kmax,1);T(1)=norm(r);k=1;if iprec==2L=ichol(A);for i=1:ncrj=L \ r;v=L' \ j;x=x+v;r=b-A*x;endj=L \ r;p=L' \j;while T(k)>E*norm(b) && kE*norm(b) && k

Script per il controllo della soluzione allo stazionario (t--->infinito):function [ unew,k,rr ] = ustaz(ndir,dir,H,f,kmax,ncr,nn,iprec )%function per il

```

calcolo della soluzione allo stazionario

imposizione delle condizioni al contorno

for i=1:ndir

j=dir(i);

H(j,j)=10^+15;

f(j)=0;

end

E=10^-8*norm(f);

chiamata del GCM

[ unew,k,T,Tvero,rvero,rr ] = gradconmod( H,f,E,kmax,ncr,nn,iprec );

semilogy(1:k,rr(1:k))

title('CONVERGENZA DEL METODO')

xlabel('ITERAZIONI')

ylabel('RESIDUO RELATIVO')

legend('MESH')

end

6) SOLUZIONE NEL TEMPO

Ora si è pronti per iniziare il ciclo temporale. La soluzione nel tempo viene calcolata con il metodo alle differenze finite, e considerando la variabile temporale come parte del termine noto. All'interno del ciclo while vengono calcolati K1,K2 e rhs del sistema lineare scaturito dall'imposizione del metodo teta per risolvere il sistema lineare nel tempo:

Prima di procedere alla risoluzione del sistema al passo temporale corrente, occorre imporre le condizioni al contorno di dirichlet che nel nostro caso sono uguali a 0 nei nodi che compongono il dominio della nostra membrana.

Per non perdere la simmetria di K1, le condizioni al contorno vengono applicate sostituendo il valore sulla diagonale della matrice con un valore penalty (ca 10^15) e l'elemento del termine noto pari a 0. Ora il sistema può essere risolto utilizzando il Gradiente Coniugato modificato essendo K1 matrice sparsa, simmetrica e definita positiva. Il vettore delle soluzione che scaturisce dal GCM lo salvo al passo temporale corrente; ed infine aggiorno le variabili. Script per il controllo della soluzione del tempo: ```html function [u]=utempo(tempo,tmax,H,P,f,teta,deltat,ndir,dir,kmax,ncr,nn,iprec,ndt) u=zeros(nn,30); while tempo<=tmax k1=(teta*H)+(P/deltat); k2=(P/deltat)-(1-teta)*H; rhs=k2*u(:,ndt)-f; for i=1:ndir j=dir(i); k1(j,j)=10^+15; rhs(j)=0; end E=10^-8*norm(rhs); [ unew,k,T,Tvero,rvero ] = gradconmod( k1,rhs,E,kmax,ncr,nn,iprec ); ndt=ndt+1; u(:,ndt)=unew; tempo=tempo+deltat; end ``` Controllo la soluzione del punto (0,0) nel

tempofigure(2)plot(1:ndt,u(11,1:ndt));title('CONTROLLO SOLUZIONE NEL TEMPO DEL PUNTO(0,0)')xlabel('n° PASSI TEMPORALI')ylabel('VALORE DI "u"')end

7) RISULTATI, ALLEGATI E CONSIDERAZIONI

Tabella di confronto fra la soluzione numerica e quella analitica per la griglia più grossolana con il calcolo dell'errore puntuale;

In questa tabella è riassunta la soluzione allo stazionario (indipendente dal tempo) per la griglia più grossolana (mesh0). Viene riportato il confronto fra la soluzione numerica e quella reale e l'errore in valore assoluto punto per punto.

Si nota che l'errore nei punti in cui la soluzione è 0, è nell'ordine di 10^-15 e quindi trascurabile. Mentre nelle soluzioni diverse da 0 l'errore è dell'ordine di 10^-2... ancora un po' troppo elevato, ma per una griglia così poco raffinata è corretto.

8) DIAGRAMMI DI CONVERGENZA IN SCALA SEMILOGARITMICA

DEL GCMPER LA SOLUZIONE DEL SISTEMA: SI PRESENTI UN DIAGRAMMA CONTUTTI I PROFILI DI CONVERGENZA RELATIVI ALLE MESH ASSEGNATE PER ILCASO STAZIONARIO OTTENUTI PARTENDO DALLA SOLUZIONE INIZIALENULLA E UTILIZZANDO SIA K^-1=F^-1(JACOBI) CHE K^-1=(L*L’)^-1(CHOLESKY). SI PONGA SULL’ASSE DELLE ASCISSE IL NUMERO DIITERAZIONI E SULL’ASSE DELLE ORDINATE LA NORMA EUCLIDEA DELRESIDUO RELATIVO;Utilizzando i due precondizionatori e tutte le mesh, si sono plottati i grafici in figura.Si nota come in entrambi i grafici al variare del tipo di mesh (griglia) di discretizzazione del dominio varia ilnumero di iterazioni necessarie per raggiungere la convergenza. Questo è dovuto al fatto che aumentandoil raffinamento della griglia aumenta anche il numero di punti da calcolare e da salvare, aumentando anchei tempi di convergenza. Così ho che la griglia più grossolana(mesh 0) è più veloce a convergere ma è menoprecisa; mentre la griglia piùraffinata(mesh4) è molto più precisa, ma con tempi di convergenza elevati. Si nota come il numero di iterazioni tra il metodo di Cholesky e quello di Jacobi cresca notevolmente. Questo è dovuto all'approssimazione della matrice di rigidezza. Con il metodo di Jacobi approssimo la matrice solo con la sua diagonale, mentre con la decomposta incompleta di Cholesky la approssimo con la sua parte triangolare bassa; per cui più simile è la matrice utilizzata nel GCM rispetto alla matrice reale, più il metodo sarà efficace e il metodo convergerà velocemente. DIAGRAMMA DI CONVERGENZA IN SCALA SEMILOGARITMICA DEL GCM PER LA SOLUZIONE DEL SISTEMA CON 1, 10, 20, 40 E 60 ITERAZIONI PRELIMINARI EFFETTUATE CON LO SCHEMA DELLE CORREZIONI RESIDUE: PER QUESTO TEST SI UTILIZZI LA MESH PIÙ FINE NEL CASO STAZIONARIO USANDO IL PRECONDIZIONATORE DI CHOLESKY. Con il metodo delle correzioni residue implementato al gradiente coniugato faccio in modo di

scegliere la soluzione iniziale del Gradiente coniugato modificato. Quanto più la soluzione sarà vicina agli autovalori del sistema tanto più rapido sarà il GCM a convergere alla soluzione reale del sistema. Notiamo come tanto più aumento il numero di correzioni residue pre GCM, tanto più il metodo convergerà velocemente. Bisogna trovare il giusto compromesso tra numero di CR e costo computazionale per rendere il metodo al massimo dell'efficienza.

TABELLA RIASSUNTIVA DEI RISULTATI DI CONVERGENZA IN STAZIONARIO, IN CUI SIANO RIPORTATI AD OGNI RAFFINAMENTO SUCCESSIVO DELLA TRIANGOLAZIONE:

(a) IL VALORE DELL'ERRORE,

(b) IL RAPPORTO FRA L'ERRORE CORRISPONDENTE ALLA GRIGLIA CORRENTE E ALLA GRIGLIA MENO RAFFINATA IMMEDIATAMENTE PRECEDENTE.

Con questa tabella si relaziona l'errore del metodo per il risultato allo stazionario. La norma euclidea dell'errore è per definizione: l'integrale può anche essere

calcolato numericamente mediante la sommatoria sugli elementi e estesa a tutti i triangoli che condividono il nodo "i": Quindi è come pesare la differenza delle soluzioni con l'area affer