[Analisi FEM]
Ciao ragazzi, mi sto per la prima volta approcciando al metodo degli elementi finiti quindi non siate crudeli.
Per prima cosa sto cercando di capire il concetto alla base del metodo:
Ogni problema reale può essere modellizzato tramite un'equazione differenziale, la cui soluzione che è una funzione come u(x), sopratutto per geometrie complesse, è difficile da ricavare. Quindi si è pensato di discretizzare il dominio di applicazione in una serie di elementi finiti, delimitati da nodi. A ciascun elemento finito si associa un sistema di equazioni algebriche, scrivibile tramite una matrice di rigidezza(Quindi un'equazione per ciascun nodo). La sovrapposizione di questi sistemi e la risoluzione alla fine mi permette di ottenere nei diversi nodi il valore della funzione u e questi valori poi vengono interpolati con un polinomio che approssima quindi la funzione u(x) da noi cercata. Se voglio una soluzione più approssimata di quella esatta dunque mi conviene aumentare gli elementi finiti e il numero di nodi. Più o meno?
Per prima cosa sto cercando di capire il concetto alla base del metodo:
Ogni problema reale può essere modellizzato tramite un'equazione differenziale, la cui soluzione che è una funzione come u(x), sopratutto per geometrie complesse, è difficile da ricavare. Quindi si è pensato di discretizzare il dominio di applicazione in una serie di elementi finiti, delimitati da nodi. A ciascun elemento finito si associa un sistema di equazioni algebriche, scrivibile tramite una matrice di rigidezza(Quindi un'equazione per ciascun nodo). La sovrapposizione di questi sistemi e la risoluzione alla fine mi permette di ottenere nei diversi nodi il valore della funzione u e questi valori poi vengono interpolati con un polinomio che approssima quindi la funzione u(x) da noi cercata. Se voglio una soluzione più approssimata di quella esatta dunque mi conviene aumentare gli elementi finiti e il numero di nodi. Più o meno?
Risposte
La sezione giusta credo sia Analisi Numerica. Sostanzialmente sì, sono particolarmente importanti quando la geometria del dominio nel quale devi risolvere il problema non è "buona" e magari la soluzione non è molto regolare. L'esempio classico è l'equazione della trave in 1D. Il principio matematico dietro, in genere, è la formulazione debole del problema, che una volta discretizzata dà origine ad un sistema lineare (se la forma è lineare).
So che a ingegneria spesso questo non viene affrontato dal questo punto di vista, ma ecco un esempio classico e introduttivo. Parti da $-u''(x)=f(x)$ con condizioni di Dirichlet omogenee (cioè $u(0)=u(1)=0$), per semplicità. Supponi $f \in L^2(0,1)=\{ f \text{ t.c. } ( \int_0^1 f^2(x) \text{dx})^2 < \infty \}$
Prendi una funzione "test" $v$, che si annulla in $x=0$ e $x=1$ in un certo spazio $V$ che tra poco specifico. Moltiplica ambo i membri per $v$, e poi integra:
$$\int_0^1 - v(x) u''(x) dx = \int_0^1 v(x)f(x)dx$$ Integra per parti il membro di sinistra, usando il fatto che $v(0)=v(1)=0$: $$\int_0^1 v'(x)u'(x)dx = \int_0^1 v(x) f(x)xdx$$
La formulazione *debole* equivale a cercare una funzione $u$ (che sta anch'essa in $V$) tale che valga
$$\int_0^1 v'(x)u'(x)dx = \int_0^1 v(x) f(x)xdx$$
per ogni $v \in V$. Lo spazio $V$ è lo spazio di Sobolev $H_0^1(0,1)$. Semplificando al massimo: una funzione che sta qui dentro può essere una **spezzata** "a tratti", proprio come immagini debba essere dal punto di vista fisico. In particolare, una funzione che sta qui dentro NON ha derivata seconda, come invece è necessario per una soluzione **forte** del problema di partenza, in quanto l'incognita $u$ aveva la derivata seconda.
Ora il punto: dopo aver partizionato l'intervallo $(0,1)$ in $N$ nodi $x_i$, lo spazio $V$ è generato dalla funzioni a cappello $phi_i(x)$ "(hat functions)", perciò siccome deve valere per ogni $v \in V$, allora deve valere per ogni elemento della base $\{ \phi_i \}_i$:
$$\int_0^1 u'(x) \phi_i'(x)dx = \int_0^1 f(x) \phi_i(x)dx, \quad \quad i=1,\ldots,N$$
Ma anche $u \in V$: perciò deve potersi scrivere come combinazione lineare degli elementi di una base, e cioè $$u(x) = \sum_j u_j \phi_j(x)$$
e dunque, inserendo questo nell'ultimo integrale:
$$ \sum_j u_j \int_0^1 \phi_j'(x) \phi_i'(x)dx = \int_0^1 f(x) \phi_i(x)dx, \quad \quad i=1,\ldots,N$$
Ma questo non è altro che la riscrittura di un sistema lineare $A \vec{u} = \vec{f}$ dove $\vec{u}$ è il vettore contenente i coefficienti $u_j$, mentre $A$ è la matrice di entrate $$( A)_{ij} = \int_0^1 \phi_j'(x) \phi_i'(x)dx $$ e invece $\vec{f}$ è un vettore di entrate $( f )_i = \int_0^1 f(x) \phi_i(x)$
Questi integrali si possono calcolare molto semplicemente, in quanto il supporto di $\phi(x)$ è $[x_{i-1},x_{i+1}]$ dando così origine ad una matrice tridiagonale. Per $f_i$ nota che $f_i = \int_0^1 f(x) \phi_i(x)dx = \int_{x_{i-1}}^{x_{i+1}} f(x) \phi(x)dx$ e questo integrale lo puoi calcolare con qualche formula di quadratura. Nota che nel caso 2D/3D la strategia è *molto* più complessa perché serve ragionare in termini di triangolazioni, e si ricorre all'utilizzo di formule di quadratura
- - - - - -
Riguardo alla tua ultima domanda: aumentare l'accuratezza si può fare in due modi:
1. Infittire la griglia, che implica un considerevole sforzo computazionale
2. Grazie ad alcune stime sull'errore, potresti usare elementi grado più alto, ad esempio quadratici,cubici,... ma ciò richiede una maggiore regolarità della soluzione. Nel momento in cui questa abbia una regolarità più bassa del grado dell'elemento, non avrai la buona approssimazione che ti aspetti.
So che a ingegneria spesso questo non viene affrontato dal questo punto di vista, ma ecco un esempio classico e introduttivo. Parti da $-u''(x)=f(x)$ con condizioni di Dirichlet omogenee (cioè $u(0)=u(1)=0$), per semplicità. Supponi $f \in L^2(0,1)=\{ f \text{ t.c. } ( \int_0^1 f^2(x) \text{dx})^2 < \infty \}$
Prendi una funzione "test" $v$, che si annulla in $x=0$ e $x=1$ in un certo spazio $V$ che tra poco specifico. Moltiplica ambo i membri per $v$, e poi integra:
$$\int_0^1 - v(x) u''(x) dx = \int_0^1 v(x)f(x)dx$$ Integra per parti il membro di sinistra, usando il fatto che $v(0)=v(1)=0$: $$\int_0^1 v'(x)u'(x)dx = \int_0^1 v(x) f(x)xdx$$
La formulazione *debole* equivale a cercare una funzione $u$ (che sta anch'essa in $V$) tale che valga
$$\int_0^1 v'(x)u'(x)dx = \int_0^1 v(x) f(x)xdx$$
per ogni $v \in V$. Lo spazio $V$ è lo spazio di Sobolev $H_0^1(0,1)$. Semplificando al massimo: una funzione che sta qui dentro può essere una **spezzata** "a tratti", proprio come immagini debba essere dal punto di vista fisico. In particolare, una funzione che sta qui dentro NON ha derivata seconda, come invece è necessario per una soluzione **forte** del problema di partenza, in quanto l'incognita $u$ aveva la derivata seconda.
Ora il punto: dopo aver partizionato l'intervallo $(0,1)$ in $N$ nodi $x_i$, lo spazio $V$ è generato dalla funzioni a cappello $phi_i(x)$ "(hat functions)", perciò siccome deve valere per ogni $v \in V$, allora deve valere per ogni elemento della base $\{ \phi_i \}_i$:
$$\int_0^1 u'(x) \phi_i'(x)dx = \int_0^1 f(x) \phi_i(x)dx, \quad \quad i=1,\ldots,N$$
Ma anche $u \in V$: perciò deve potersi scrivere come combinazione lineare degli elementi di una base, e cioè $$u(x) = \sum_j u_j \phi_j(x)$$
e dunque, inserendo questo nell'ultimo integrale:
$$ \sum_j u_j \int_0^1 \phi_j'(x) \phi_i'(x)dx = \int_0^1 f(x) \phi_i(x)dx, \quad \quad i=1,\ldots,N$$
Ma questo non è altro che la riscrittura di un sistema lineare $A \vec{u} = \vec{f}$ dove $\vec{u}$ è il vettore contenente i coefficienti $u_j$, mentre $A$ è la matrice di entrate $$( A)_{ij} = \int_0^1 \phi_j'(x) \phi_i'(x)dx $$ e invece $\vec{f}$ è un vettore di entrate $( f )_i = \int_0^1 f(x) \phi_i(x)$
Questi integrali si possono calcolare molto semplicemente, in quanto il supporto di $\phi(x)$ è $[x_{i-1},x_{i+1}]$ dando così origine ad una matrice tridiagonale. Per $f_i$ nota che $f_i = \int_0^1 f(x) \phi_i(x)dx = \int_{x_{i-1}}^{x_{i+1}} f(x) \phi(x)dx$ e questo integrale lo puoi calcolare con qualche formula di quadratura. Nota che nel caso 2D/3D la strategia è *molto* più complessa perché serve ragionare in termini di triangolazioni, e si ricorre all'utilizzo di formule di quadratura
- - - - - -
Riguardo alla tua ultima domanda: aumentare l'accuratezza si può fare in due modi:
1. Infittire la griglia, che implica un considerevole sforzo computazionale
2. Grazie ad alcune stime sull'errore, potresti usare elementi grado più alto, ad esempio quadratici,cubici,... ma ciò richiede una maggiore regolarità della soluzione. Nel momento in cui questa abbia una regolarità più bassa del grado dell'elemento, non avrai la buona approssimazione che ti aspetti.
Dimenticavo, una volta calcolato $\vec{u}$, soluzione del sistema lineare precedente, la soluzione è data da $$u(x)= \sum_j u_j \phi_j(x)$$, ossia l'interpolante lineare a tratti. Nel tuo post parli della somma dei contributi locali, questo è un passaggio fondamentale perché in generale la matrice non è semplice come quella sopra, ad esempio se il problema è in 2D, e di fatto la si "assembla" sommando dei contributi locali. Per i dettagli ti consiglio il libro di Quarteroni "Modellistica Numerica per problemi differenziali" ,oppure per un'introduzione più ingegneristica dove c'è anche un po' di codice MatLab: The Finite Element Method: Theory, Implementation, and Applications (qui il pdf: http://matematicas.unex.es/~coco/Modelizacion/FEMLarson_Bengzon.pdf), vedi capitolo 4 dove fanno il problema di Poisson in 2D
Accedi a tutti gli appunti
Tutor AI: studia meglio e in meno tempo