Un paio d`esempi completi con FreeFem++ Un paio d`esempi

Un
Un paio
paio d'esempi
d'esempi completi
completi con
con
FreeFem++
FreeFem++
Isacco Perin
University of Verona
Master's Degree in Mathematics and Applications
Isacco Perin (University of Verona)
13/05/2013
FreeFem++
Introduzione
L'obiettivo di tale presentazione è di fornire, mediante un paio
d'esempi, le linee guida per programmare in FreeFem++.
FreeFem++
FreeFem++ èè un
un codice
codice per
per la
la risoluzione
risoluzione di
di PDE
PDE mediante
mediante ilil
metodo
metodo degli
degli elementi
elementi finiti.
finiti.
Come linguaggio di programmazione, esso ha una sintassi di base
simile al C++, pertanto è possibile dichiarare variabili e funzioni
all'interno di uno script FreeFem++, eseguire operazioni aritmetiche,
utilizzare cicli o istruzioni condizionali.
Isacco Perin (University of Verona)
13/05/2013
FreeFem++
Introduzione
Codice e software sono stati sviluppati da
O. Pìronneau, F. Hecht, e A. Le Hyaric.
Il software e il manuale d'uso sono scaricabili gratuitamente da
internet al link http://www.freefem.org
Isacco Perin (University of Verona)
13/05/2013
Membrana Elastica
Definizione del Problema
Mediante la definizione di opportune PDE è possibile esprimere in
termini matematici situazioni fisiche. Dunque, al fine di mostrare
l'applicabilità di un algoritmo scritto in FreeFem++ alla realtà, gli
esempi trattati in seguito si baseranno su di una specifica situazione
fisica.
Problema:
Problema: Membrana
Membrana Elastica
Elastica
Come reagisce una membrana elastica attaccata ad un supporto
rigido se sottoposta ad un carico?
Dal punto di vista matematico possiamo esprimere tale situazione
nel modo seguente:
Isacco Perin (University of Verona)
13/05/2013
Definizione dal punto di vista Matematico
Definizione del Problema
Consideriamo una membrana elastica Ω attaccata ad un supporto
rigido Γ , dunque si ha che ∂ Ω=Γ .
Consideriamo una forza, f ( x) dx ,esercitata su ogni elemento della
superficie dx=dx 1 dx 2 .
Per ottenere lo spostamento verticale della membrana , φ(x) ,si deve
risolvere la seguente equazione di Laplace
−Δ φ( x)= f ( x) , x ∈Ω.
Isacco Perin (University of Verona)
13/05/2013
Condizioni al Bordo
Definizione del Problema
Condizioni
Condizioni al
al bordo
bordo
Per quanto riguarda le condizioni da imporre al bordo, esse dipendono
dalla situazione fisica trattata. Di seguito presentiamo delle situazioni
in cui possiamo adottare condizioni di Dirichlet o di Neumann
omogenee.
-
La membrana è fissata al supporto planare. Allora, la condizione al
bordo da adottare è di Dirichlet omogenea ed è la seguente:
φ( x)=0, x∈Γ .
-
Il supporto non è planare, ma è un elevazione che si esprime
formalmente con z( x). Allora, la condizione al bordo da adottare è
di Dirichlet non omogenea ed è la seguente:
φ( x)= z( x) , x ∈Γ .
Isacco Perin (University of Verona)
13/05/2013
Condizioni al Bordo
Definizione del Problema
-
La membrana è fissata ad una parte del supporto Γ1 mentre è
lasciata semplicemente appesa nella parte restante, Γ 2 . Allora, per
la rigidità della membrana, le condizioni al bordo da adottare sono
di Neumann omogenee lungo Γ 2 e di Dirichlet omogenee lungo Γ1 .
φ( x)=0, x∈Γ 1 ,
∂ φ( x)
=0, x∈Γ2 .
∂n
Ovviamente, in tal caso, il bordo è Γ=Γ1∪Γ 2 .
Isacco Perin (University of Verona)
13/05/2013
Esistenza e Unicità
Problema di Laplace
Esistenza
Esistenza ee Unicità
Unicità
L'equazione di Laplace precedentemente definita, con opportune
ipotesi, ammette un'unica soluzione.
ATTENZIONE:
Ricordiamo che FreeFem++ adotta il metodo degli elementi finiti ed
esso, per qualsiasi tipo di elementi finiti si scelgano, è basato sulla
formulazione debole di una PDE.
Infatti, FreeFem++ risolve equazioni alle derivate parziali espresse in
formulazione debole.
Dunque, prima di mostrare l'algoritmo adottato per risolvere il problema
della membrana elastica, dobbiamo mostrare la formulazione debole di
tale problema.
Isacco Perin (University of Verona)
13/05/2013
Formulazione Debole
Problema di Laplace
La formulazione forte (o classica) del problema della membrana
elastica con condizioni al bordo di Dirichlet omogenee è
−Δ φ( x)= f ( x) , x ∈Ω
φ(x)=0, x ∈Γ .
La formulazione debole di tale PDE è la seguente
∫Ω ∇ φ ∇ v=∫Ω f v , ∀ v∈H (Ω) ,
1
dove H 1 (Ω)= { v∈ L2 (Ω) , ∇ v∈( L2 (Ω))×(L 2 (Ω)) } .
Isacco Perin (University of Verona)
13/05/2013
Esempio
Equazione di Laplace con condizioni di Dirichlet omogenee
Presentiamo, dunque, il primo dei due esempi.
Esempio
Esempio
Determinare la soluzione dell'equazione di Laplace in FreeFem++ nel
caso in cui Ω è un ellisse avente lunghezza del semiasse maggiore 2 e
di quello minore 1, f (x)=1 e con condizioni al bordo di Dirichlet
omogenee.
L'algoritmo è riportato di seguito, con opportuni commenti.
Isacco Perin (University of Verona)
13/05/2013
Esempio
Equazione di Laplace con condizioni di Dirichlet omogenee
// Declaration of some global parameters
real theta=4.*pi/3., a=2., b=1.;
// Force function
func f = 1.0;
// Defining the boundary as an ellipse
border gamma1(t=0,2*pi) {x=a*cos(t); y=b*sin(t);}
// Mesh
mesh M = buildmesh(gamma1(80));
plot(M, wait=true, ps="Mesh1.eps");
Osservazione: con il comando 'buildmesh' viene costruita una
triangolazione del dominio.
Isacco Perin (University of Verona)
13/05/2013
Esempio
Equazione di Laplace con condizioni di Dirichlet omogenee
// Build the finite element space using P2 elements
fespace Vh(M,P2);
// Define u and w as piecewise-P2 function
Vh uh, w;
// Variational Formulation of our equation
problem Laplace(uh,w) =
int2d(M)(dx(uh)*dx(w)+dy(uh)*dy(w))
- int2d(M)(f*w)
+ on(gamma1,uh=0);
// Dirichlet boundary condition
Osservazione: la PDE viene discretizzata adottando il metodo degli
elementi finiti del secondo ordine sulla triangolazione del dominio.
// Solve the PDE
Laplace;
Isacco Perin (University of Verona)
13/05/2013
Esempio
Equazione di Laplace con condizioni di Dirichlet omogenee
// Plot of the solution
plot(uh, wait=true, ps="Solution1.eps");
Isacco Perin (University of Verona)
13/05/2013
Grafici
Equazione di Laplace con condizioni di Dirichlet omogenee
Di seguito sono riportati i grafici ottenuti.
La figura a sinistra rappresenta la mesh adottata, mentre quella a destra
rappresenta le curve di livello relative alla soluzione numerica della PDE
data.
Isacco Perin (University of Verona)
13/05/2013
Procedura generale
Programmare in FreeFem++
La procedura seguita in tale algoritmo, può essere generalizzata nel
modo seguente:
-
Definisco il bordo del dominio;
-
Costruisco la mesh;
-
Definisco lo spazio su cui applico il metodo degli elementi finiti;
-
Definisco le opportune funzioni su tale spazio;
-
Inserisco la formulazione variazionale della PDE che voglio risolvere.
Isacco Perin (University of Verona)
13/05/2013
Esempio
Verifica dei risultati ottenuti nel primo esempio
Esempio
Esempio
Nel secondo esempio verifichiamo i risultati ottenuti nel primo.
L'idea è di modificare la situazione precedente in modo tale da
considerare una PDE di cui conosciamo la soluzione esatta per poterla
confrontare con quella che si ottiene con FreeFem++.
Consideriamo pertanto come bordo il circolo unitario e come forzante la
funzione
2
2
2
2
2
2
f ( x)=−4(cos(x + y −1)−( x + y )sin (x + y −1)).
Isacco Perin (University of Verona)
13/05/2013
Esempio
Verifica dei risultati ottenuti nel primo esempio
Come condizioni al bordo consideriamo delle condizioni miste: di
Neumann non omogenee lungo una parte del bordo, mentre di Dirichlet
omogenee lungo la rimanente parte.
Sotto le ipotesi fatte, la soluzione esatta dell'equazione di Laplace è data
da
2
2
φ( x)=sin (x + y −1).
2
Dunque, nel seguente algoritmo calcoleremo l'errore in norma L tra la
soluzione esatta e quella ottenuta in FreeFem++,
ε=∥φ−φF∥L .
2
Isacco Perin (University of Verona)
13/05/2013
Esempio
Verifica dei risultati ottenuti nel primo esempio
Inoltre, stimeremo l'ordine di convergenza del metodo calcolando prima
l'errore ottenuto con n punti di discretizzazione lungo il bordo e poi
quello ottenuto con 2 n punti.
La stima si ottiene calcolando il logaritmo del rapporto di tali errori.
Di seguito è riportato l'algoritmo scritto in FreeFem++.
Isacco Perin (University of Verona)
13/05/2013
Esempio
Verifica dei risultati ottenuti nel primo esempio
// Declaration of some global parameters
real theta=4.*pi/3., a=1., b=1.;
real[int] L2error(2); // an array with 2 values
// Force function
func f=-4*(cos(x^2+y^2-1)-(x^2+y^2)*sin(x^2+y^2-1));
// Exact Solution
func u=sin(x^2+y^2-1);
// Defining the boundary as an ellipse
border gamma1(t=0,theta) {x=a*cos(t); y=b*sin(t);}
border gamma2(t=theta,2*pi) {x=a*cos(t); y=b*sin(t);}
plot(gamma1(50)+gamma2(50), wait=true, ps="Border.eps");
Isacco Perin (University of Verona)
13/05/2013
Esempio
Verifica dei risultati ottenuti nel primo esempio
// Evaluate the error between the exact solution and the
solution calculated with FreeFem++
for(int n=0;n<2;n++)
{
// Mesh
mesh M=buildmesh(gamma1(20*(n+1))+gamma2(10*(n+1)));
// Build the finite element space using P2 elements
fespace Vh(M,P2);
// Define u and phi as piecewise-P2 function
Vh uh, w;
Isacco Perin (University of Verona)
13/05/2013
Esempio
Verifica dei risultati ottenuti nel primo esempio
// Variational Formulation of our equation
problem Laplace(uh,w) =
int2d(M)(dx(uh)*dx(w)+dy(uh)*dy(w))
- int2d(M)(f*w)
- int1d(M,gamma2)(2*w) // Neumann boundary condition
+ on(gamma1,uh=0);
// Dirichlet boundary condition
// Solve the PDE
Laplace;
plot(M, wait=true, ps="Mesh.eps");
plot(uh, wait=true, ps="NumericalSolution.eps");
L2error[n]= sqrt(int2d(M)((uh-u)^2));
}
Isacco Perin (University of Verona)
13/05/2013
Esempio
Verifica dei risultati ottenuti nel primo esempio
// Show the error and the rate convergence
for(int n=0;n<2;n++)
cout << " L2error " << n << " = "<< L2error[n] << endl;
cout << " convergence rate = "<<
log(L2error[0]/L2error[1])/log(2.) << endl;
Isacco Perin (University of Verona)
13/05/2013
Grafici
Verifica dei risultati ottenuti nel primo esempio
Di seguito sono riportati i grafici ottenuti.
La figura a sinistra rappresenta la mesh adottata, mentre quella a destra
rappresenta le curve di livello relative alla soluzione numerica della PDE
data.
Isacco Perin (University of Verona)
13/05/2013
Risultati
Verifica dei risultati ottenuti nel primo esempio
In output si ottengono i seguenti risultati:
L2error 0 = 0.018358
L2error 1 = 0.0046485
convergence rate = 1.98157
Isacco Perin (University of Verona)
13/05/2013
Bibliography
[1] F.Hecht , FreeFem++
[2] O. Pìronneau, FreeFEM User Manual (2001)
Isacco Perin (University of Verona)
13/05/2013