Relazione esercizio 2 Esame: Metodi Numerici

U=L';

y_chol=L\l;

sol_chol=U\y_chol;

end

%controllo sulle condizioni al contorno

if (bc_l_t == 'D')

sol_chol = [bc_l_v;sol_chol];

end

if (bc_r_t == 'D')

sol_chol = [sol_chol;bc_r_v];

end

t_chol=toc;

Tra i metodi diretti il migliore è sicuramente l’algoritmo di Cholesky in quanto è un

algoritmo stabile, utilizza un solo elemento per la fattorizzazione e riduce della metà la

complessità computazionale rispetto agli altri metodi (ad esempio Gauss). Il metodo di

Cholesky però risulta valido per matrici simmetriche e definite positive. Difatti, durante

le lezioni si è affermato spesso che Cholesky oltre ad essere un solutore è anche un

verificatore per le matrici in quanto, se non funziona, la matrice non soddisfa le

condizioni di simmetria o di positività. In ogni caso MatLab restituirà un messaggio di

errore nella Command Window. In tutte le mesh realizzate la matrice di rigidezza è

risultata sempre simmetrica e definita positiva. Per rispondere al quesito l’algoritmo di

Cholesky è stato confrontato sulla base del tempo con tutti i metodi iterativi inseriti nel

codice, grazie al comando tic toc di MatLab. Partendo dallo stesso valore della

tolleranza e dallo stesso numero massimo di iterazioni, per ogni mesh i risultati ottenuti

sono quelli riportate nelle tabelle seguenti: 11

Relazione Progetto –Corso di Metodi Numerici per l’Ingegneria Civile Magistrale- A.A. 2016-2017 – Esercizio 2

iterm_max= 940

toll=10^-7 12

Relazione Progetto –Corso di Metodi Numerici per l’Ingegneria Civile Magistrale- A.A. 2016-2017 – Esercizio 2

Come si può notare i metodi di Jacobi e Gauss-Seidel non convergono entro il numero

massimo di iterazioni imposto all’aumentare della discretizzazione della mesh. Mentre

i metodi SOR e del Gradiente coniugato, pur andando a convergenza rapidamente,

impiegano sicuramente un tempo più alto rispetto al metodo diretto di Cholesky. 13

Relazione Progetto –Corso di Metodi Numerici per l’Ingegneria Civile Magistrale- A.A. 2016-2017 – Esercizio 2

Punto 5

Tenendo presente che l’equazione differenziale in (1) ammette come integrale

generale la funzione:

dove C0 e C1 sono delle costanti da determinarsi imponendo le condizioni al bordo,

grafico la soluzione esatta e quella approssimata ottenuta per

si riportino in un

successivi raffinamenti della mesh.

La forma della soluzione esatta del problema complesso è stata inserita nella Functions

fornita per il problema semplificato ue_b modificandola per farla funzionare nello

script inerente il problema (1). Successivamente a mano si sono determinate le costanti

C0 e C1 sostituendo i valori delle condizioni di Dirichlet fornite sul testo. Le costanti

valgono rispettivamente: C0=-4.215953442;

C1=-5.784046558;

All’interno dello script questo punto è stato implementato con la seguente metodologia:

figure(4);%grafico soluzione esatta e soluzione approssimata

fplot(@(x) ue(x),[0 1],'r');

x_plot=linspace(0,1,length(uh));

plot(x_plot,uh,'--','LineWidth',2);

title('soluzione esatta e soluzioni approssimate metodo Cholesky');

xlabel('lunghezza fune');

ylabel('tensione');

legend('soluzioni approssimate','soluzione esatta');

Si può osservare l’utilizzo del comando fplot, in quanto avendo già l’espressione della

soluzione esatta attraverso l’apposita Function basta richiamare la funzione e inserire il

dominio nel quale la si vuole plottare. Il grafico risultante dalle seguenti operazioni è il

seguente: 14

Relazione Progetto –Corso di Metodi Numerici per l’Ingegneria Civile Magistrale- A.A. 2016-2017 – Esercizio 2

Punto 6

Per successivi raffinamenti della mesh, si ricavi la derivata della soluzione calcolata

e la si confronti graficamente con la derivata della soluzione esatta.

Per ricavare la derivata della soluzione approssimata, è stata utilizzata la formula delle

differenze in avanti (limite del rapporto incrementale). Per il calcolo quindi, con le

formule delle differenze in avanti, sono necessari 2 punti: il punto nel quale si vuole

valutare e il punto successivo a quello in esame e nel caso dell’esercizio sono entrambi

noti. La distanza tra due nodi successivi è stata scelta come incremento h del rapporto

incrementale. Si riporta parte di codice dove si esegue quanto affermato:

figure(3);%grafico derivata esatta e derivata approssimata

fplot(@(x) ue_prime(x),[0,1],'k','LineWidth',2);%derivata esatta

for i=1:length(sol_chol)-1

ue_prime_app(i)=(sol_chol(i+1)-sol_chol(i))./h; );%differenze avanti

end

%grafico derivata approssimata

x_plot=linspace(0,1,length(ue_prime_app));

plot(x_plot,ue_prime_app,'--g');

title( 'derivata esatta e derivate approssimate');

xlabel('derivate');

ylabel('lunghezza fune');

legend('derivate approssimate','derivata esatta');

hold on

grid on 15

Relazione Progetto –Corso di Metodi Numerici per l’Ingegneria Civile Magistrale- A.A. 2016-2017 – Esercizio 2

Punto 7

Per ciascuna delle mesh considerate, si calcoli l’errore commesso rispetto alla

soluzione esatta e rispetto alla sua derivata. In entrambi i casi si realizzi un grafico

dell’errore, utilizzando una opportuna scala per la rappresentazione.

È stato richiesto di calcolare l’errore tra la soluzione approssimata e la soluzione esatta,

pertanto si è dovuta definire la soluzione esatta su ogni punto del dominio in cui è

presente un nodo nella discretizzazione (per ogni raffittimento della Mesh).

Successivamente per valutare l’errore assoluto tra le due soluzioni si è calcolata la

norma della differenza tra i valori della soluzione approssimata e i valori della soluzione

esatta. È stata scelta la norma Euclidea, già presente in matlab con il comando norm,

che calcola rispettivamente la radice quadrata della somma del quadrato delle

componenti di un vettore in modo da avere così una visione dell’errore “totale” per ogni

mesh. In generale l’utilizzo dell’errore assoluto rispetto a quello relativo è giustificato

dal fatto che nel seguente problema conosciamo la soluzione esatta. Di seguito sono

visualizzabili i risultati nel grafico seguente che mostra in ascissa i nodi della mesh

considerata e in ordinata (in scala logaritmica grazie al comando semilogy) i valori

dell’errore assoluto rispetto alla soluzione esatta.

figure(5);%grafico errore assoluto soluzione esatta e soluzioni approssimate

%x=linspace(0,1,Nh);

ue_esatta=ue_b(P');

err_ass_funzione=norm(sol_chol-ue_esatta);

semilogy(Nh,err_ass_funzione,'r*','MarkerSize',9);

title('errore soluzione');

xlabel('raffittimento mesh'); 16

Relazione Progetto –Corso di Metodi Numerici per l’Ingegneria Civile Magistrale- A.A. 2016-2017 – Esercizio 2

ylabel('errore assoluto');

legend('errore assoluto');

hold on

grid on

Si può osservare come all’aumentare della discretizzazione l’errore rispetto alla

soluzione esatta diminuisca (quello che ci aspettavamo dato che stiamo utilizzando un

approccio agli elementi finiti). Per calcolare invece l’errore tra la derivata della

soluzione approssimata e la derivata della soluzione esatta, la derivata della soluzione

esatta è stata calcolata su ogni punto del dominio in cui è presente un nodo (per ogni

raffittimento della Mesh). Anche in questo caso è stata scelta la norma Euclidea. Di

seguito i risultati sono stati riportati in un grafico che mostra in ascissa i nodi della mesh

considerata e in ordinata (in scala logaritmica grazie al comando semilogy) i valori

dell’errore assoluto rispetto alla derivata della soluzione esatta.

figure(6);% errore assoluto derivata esatta e derivate approssimate

ue_prime_esatta=ue_prime_b(P,bc_l_t,bc_r_t,bc_l_v,bc_r_v);

err_ass_derivata=norm(ue_prime_app-ue_prime_esatta);

semilogy(Nh,err_ass_derivata,'r*','MarkerSize',5); 17

Relazione Progetto –Corso di Metodi Numerici per l’Ingegneria Civile Magistrale- A.A. 2016-2017 – Esercizio 2

Si nota che in questo caso all’aumentare del numero dei nodi l’errore rispetto alla

derivata esatta diminuisce. Nel caso del caso della derivata gli errori sono più grandi

rispetto a quelli ottenuti per le soluzioni. Questa differenza è dovuta alla formula alle

differenze in avanti per il calcolo della derivata approssimata che restituisce valori poco

precisi.

Punto 8

Si modifichi lo script in riferimento al calcolo del vettore dei carichi. In particolare

si utilizzi il metodo numerico più appropriato tra quelli visti per il calcolo degli

integrali che costituiscono gli elementi di tale vettore motivando nella relazione la

scelta effettuata.

Tutte le operazioni effettuate fino ad ed adesso sono state calcolate con un vettore dei

carichi ricavato attraverso il metodo di Simpson-Cavalieri che troviamo sul codice già

implementato. Come sappiamo dalla teoria la formula di Simpson-Cavalieri è un caso

particolare delle formule di Newton-Cotes che altro non sono che delle formule di

quadratura interpolatorie su nodi equi spaziati. Si è scelto di non stravolgere troppo il

codice e di voler calcolare il vettore dei carichi implementando il metodo del 3/8 dopo

aver fatto delle prove sulle formule di quadratura interpolatorie. Siccome troveremo più

punti ci aspetteremo di ottenere un vettore dei carichi più preciso con un errore minore.

Pertanto, per visualizzare questi aspetti e osservare cosa accade, il codice è stato

modificato nel seguente modo: 18

Relazione Progetto –Corso di Metodi Numerici per l’Ingegneria Civile Magistrale- A.A. 2016-2017 – Esercizio 2

%calcolo dei punti delle funzioni interpolanti

el_qns2 =

[el_ns_crd(1),(0.33*(el_ns_crd(2)el_ns_crd(1))+el_ns_crd(1)),(0.66*(el_ns_cr

d(2)-el_ns_crd(1))+el_ns_crd(1)),el_ns_crd(2)];

for the local load vector

el_phi1s = P1_phi1(el_qns2,el_ns_crd(1),el_ns_crd(2));

el_phi2s = P1_phi2(el_qns2,el_ns_crd(1),el_ns_crd(2));

% evaluate right-hand side function

el_fs = fr(el_qns2);

% calcolo vettore locale dei carichi

l_k = zeros(2,1);%inizializzo il vettore dei carichi locale

el_F1s = el_fs .* el_phi1s

el_F2s = el_fs .* el_phi2s

l_k(1) = el_h * (el_F1s(1) + 3*el_F1s(2) + 3*el_F1s(3)+el_F1s(4))/8;

l_k(2) = el_h * (el_F2s(1) + 3*el_F2s(2) + 3*el_F2s(3)+el_F2s(4))/8;

I risultati ottenuti hanno evidenziato il seguente dato: se si fa eseguire il codice con la

modifica apportata al vettore dei carichi si può visualizzare il grafico dell’errore rispetto

alla soluzion