Problema di Dirichlet

TRIANGOLARE INFERIORE

SOSTITUZIONE IN AVANTI

Nel caso di sistemi Ax=b non triangolari, si cerca di ricondursi a più sistemi triangolari equivalenti

fattorizzando la matrice A.

• Fattorizzazione A=LU

tipo Doolittle:

Utilizzando il metodo di eliminazione di Gauss:

iterando il procedimento per j=1..n si sottrae a ogni riga i (con i=j+1,..n) la riga j moltiplicata per un

coefficiente mij. Nel caso della prima riga (j=1) si ha:

Se un elemento sulla diagonale è nullo il procedimento si arresta. (LU no pivoting)

Pagina 10 di 29

In questo caso occorre permutare righe o

colonne della matrice A in modo da non avere

elementi diagonali nulli. Si ottiene la

fattorizzazione PAQ=LU con P e Q matrici di

permutazione. (LU pivoting totale)

Tuttavia questo metodo è troppo oneroso dal

punto di vista computazionale e si preferisce

utilizzare il metodo del pivoting parziale che,

solo quando incontra il termine diagonale

nullo, permuta righe o colonne della parte di

matrice non eliminata.

PA=LU (LU pivoting parziale)

3

Costo computazionale n /3

T

• Fattorizzazione di Cholesky A=LL

+

Valida solo per A Sym 3

Costo computazionale n /6

Si riporta la funcion utilizzata per testare i diversi metodi: Pagina 11 di 29

Dai risultati ottenuti nella risoluzione del sistema lineare al variare della mesh si può notare che

l’errore è nullo o comunque vicino alla precisione di macchina e aumenta con il numero di

operazione effettuate. I tempi, dipendenti dal numero di operazioni effettuate, aumentano

all’aumentare della dimensione del sistema e la fattorizzazione ha sempre un costo

computazionale maggiore della risoluzione dei sistemi. Pagina 12 di 29

All’aumentare della dimensione la differenza tra i metodi è più marcata e permette di verificare che

Cholesky, quando applicabile, è il miglior metodo diretto di risoluzione. Tuttavia la fattorizzazione

per matrici sparse come nel nostro caso può portare al fenomeno del fill in, ovvero a un’aumento

dei termini non nulli della matrice. Pagina 13 di 29

- Metodi iterativi

Permettono di ottenere, sotto opportune ipotesi, la convergenza alla soluzione esatta in un numero

infinito di passi. Dovendo interrompere il procedimento a un numero massimo di iterazioni

prefissato si avrà quindi, oltre agli errori di macchina, una soluzione approssimata. Questi metodi

lasciando inalterata la struttura della matrice permettono di sfruttarne la sparsità e diventano

particolarmente vantaggiosi all’aumentare delle dimensioni del sistema.

Decomponendo A in modo da

ottenere poi un sistema di

facile risoluzione, ad esempio

con M diagonale non singolare

partendo da un vettore iniziale x0 si ottiene

una successione xk 2

con costo computazionale di circa n a ogni iterazione.

considerando

l’errore al passo k

si ha convergenza se il raggio spettrale, autovalore di modulo più grande, vale

che indicando la riduzione dell’errore mostra anche quanto velocemente converge il metodo.

Tuttavia a volte è oneroso da calcolare e si preferisce utilizzare cautelativamente la norma naturale

di B che risulta essere sempre maggiore del raggio spettrale e più semplice da determinare. Può

comunque capitare per matrici fortemente malcondizionate come nel nostro caso che il metodo

non converga. A seconda della decomposizione della matrice A e della conseguente matrice di

iterazione B si hanno i metodi: • Jacobi

spostamenti simultanei

• Gauss Seidel

spostamenti successivi

Dove con il metodo di Gauss Seidel a ogni iterazione nel calcolo di xi si utilizzano le componenti xj

già calcolate nella stessa iterazione; in questo modo, se si ha la convergenza, la velocità è

maggiore rispetto il metodo di Jacobi. Pagina 14 di 29

• Metodo SOR: utilizzato per ottenere la convergenza, o velocizzare se convergente:

Dove si ha convergenza se il raggio spettrale di H è

minore di 1 e in particolare si può velocizzare

determinando il valore di omega lo rende minimo:

Metodo del Gradiente Coniugato per matrici simmetriche e definite positive:

•

si cerca di trovare il punto di minimo y della forma

quadratica il cui gradiente è equivalente al nostro

sistema:

generando una successione di valori

dove:

Utilizzando come criterio d’arresto:

- 1000 iterazioni

- 10^-6 di tolleranza sul residuo

Si riportano i risultati ottenuti risolvendo il sistema al variare della mesh: Pagina 15 di 29

Pagina 16 di 29

In tutti i casi il metodo di Gauss Seidel è più veloce del metodo di Jacobi, mentre il SOR con

omega ottimale che accelera la convergenza è migliore di entrambi; sia in termini di velocità che di

residuo.

Per 256 elementi nessuno dei metodi converge in 1000

iterazioni ed entro la tolleranza scelta eccetto il Gradiente

Coniugato che dimostra essere il miglio metodo iterativo e

con risultati paragonabili a Cholesky per i metodi diretti.

Risultati Pagina 17 di 29

- Soluzione approssimata

Si ottiene sempre la soluzione esatta nei nodi tramite la risoluzione del sistema, mentre

l’andamento approssimato sull'elemento e definito dalla funzione di forma lineare. Si può vedere

che all’aumentare del numero degli elementi si ha la convergenza con la soluzione esatta.

Pagina 18 di 29

- Derivata approssimata

Si può notare che anche per la derivata si ha la convergenza, ma questa a differenza della

soluzione non è in generale esatta nei nodi. Pagina 19 di 29

Essendo la soluzione nota solo in un numero limitato di punti, per determinare la derivata si è

utilizzato il metodo alle differenze finite. In particolare centrate per i punti intermedi e in avanti

(indietro) per il primo (ultimo) nodo utilizzando le formule di ordine superiore per migliorare la

precisione.

Dei risultati ottenuti utilizzando le differenze in avanti/indietro del primo ordine si può vedere come

l’errore sia maggiore in corrispondenza dei nodi estremi: Pagina 20 di 29

- Errori

Interpolando linearmente soluzione e derivata note solamente nei nodi è possibile calcolare l’errore

Dai grafici si può vedere che soluzione e derivata convergono molto velocemente e quest’ultima

con un’errore maggiore per via del metodo di approssimazione alle differenze finite.

Pagina 21 di 29

Integrale vettore dei carichi

Con il metodo agli elementi finiti per determinare le componenti del sistema, matrice di rigidezza K

e vettore dei carichi appaiono diversi integrali. Per quanto riguarda K vengono coinvolte le funzioni

di forma N e le relative derivate; e considerando i dati del nostro problema ed elementi lineari, le

componenti da integrare sono delle costanti per il problema semplificato e dei polinomi di secondo

grado per i termini aggiuntivi del problema generale. Come detto in precedenza il codice utilizza la

formula di quadratura di Cavalieri Simpson per calcolare l’integrale approssimato come:

Che fà parte del più generale metodo di Newton Cotes, che approssima la funzione con un

polinomio di grado n-1 e calcola l’integrale come sommatoria della funzione valutata in n nodi

equispaziati dopo averla moltiplicata per opportuni pesi. La formula di Cavalieri Simpson in

particolare approssima la funzione, utilizzando 3 punti, con un polinomio di secondo grado e quindi

è in grado di fornire l’integrale esatto per le componenti della matrice K.

Per quanto riguarda il vettore dei carichi l’integrale coinvolge la funzione di forma N lineare

moltiplicata per il carico w che in questo caso è una sinusoide, ma in generale può variare essendo

un dato del problema. Per generalizzare il codice in modo da renderlo affidabile anche per

differenti dati iniziali si è modificato lo script in riferimento al calcolo del vettore dei carichi:

Dove si è scelto la formula di quadratura di Gaussiana composita, che considera punti non

equispaziati per evitare i problemi di instabilità che si verificano con Newton Cotes all’aumentare

del numero di nodi e con n punti calcola l’integrale esatto con un ordine di 2n-1.

Visto che in generale l’errore di quadratura si riduce al ridursi

dell’intervallo, che coincide con la dimensione dell’elemento, e per via

dell’elevato grado di precisione della formula di Gauss; si è scelto di

utilizzare un unico intervallo e 4 nodi modificabili nei dati iniziali Pagina 22 di 29

Differenti dati iniziali

Aumentando l’elasticità k=100 l’abbassamento ha un’andamento simile al carico w(x)

Mentre incrementando la tensione T=10 si torna ad un risultato simile hai precedenti con

abbassamenti molto ridotti. Pagina 23 di 29

Cedimento all’estremo di destra, u(b)=-3

Aumentando la lunghezza della fune l’andamento a sinusoide del carico diventa trascurabile e

assimilabile a un carico costante w=0.5, si ottiene una soluzione simmetrica. Pagina 24 di 29

Listato Matlab

%————————————————————————————————————%

% (1) Problema generale -u''(x) +(k/T)u(x)= w(x)/T a < x < b %

% %

% u(a) = bcL =0 u(b) = bcR = 0 w(x)=1+sin(4*pi*x) a=0 b=1 %

%————————————————————————————————————%

close all;

clear all;

clc;

format short;

%1) PROBLEM DEFINITION

% domain definition

a = 0;

b = 1;

x=linspace(a,b,100);

% parametri

T=1;

k=0.1;

% termine noto w(x)

fr = @(x)carico_wT(x,T);

% boundary conditions (left and right)

% types Dirichlet, values:

bcL = 0;

bcR = 0;

%costanti da BC per u_esatta

[ c0,c1 ] = BC_1( a,b,bcL,bcR,T,k );

% exact solution and exact solution derivative:

ue = @(x)ue1( x,T,k,c0,c1 );

ue_prime = @(x)ue1_prime( x,T,k,c0,c1 );

% number of meshes

m_n = 6; %subplot pari

% quadratura

int_comp=1;

pt_gauss=4;

%------------------------------------------------------------------------%

% FEM FOR THE NUMERICAL APPROXIMATION OF THE

% SOLUTION OF THE FOLLOWING 1-D ELLIPTIC

% BVP (BOUNDARY VALUE PROBLEM)...POISSON'S EQUATION:

% NUMERICAL APPROXIMATION BY FEM:

% we compute the numerical solution by using uniform linear finite

% elements, that is we use a uniform mesh (elements of uniform size)

% and a finite-dimensional function space (where to search for the

% solution) given by continuous piecewise linear functions.

%

% To study experimentally the convergence of the FEM (that is, the numerical