ANALISI NUMERICA
Metodo di Galerkin ed introduzione agli elementi finiti
a.a. 2012–2013
Maria Lucia Sampoli
ANALISI NUMERICA – p.1/22
METODO di GALERKIN
Tutti i problemi visti possono essere ricondotti alla formulazione debole:
trovare u ∈ V
t.c. a(u, v) = F (v),
∀v∈V
con V ⊂ H 1 (Ω), spazio di Hilbert, a : V × V → R forma bilineare continua e
coerciva, F : V → R forma lineare limitata.
Formulazione approssimata: Sia Vh ⊂ V una famiglia di sottospazi (al variare di h) a
dimensione finita dim Vh = Nh :
trovare uh ∈ Vh t. c. a(uh , vh ) = F (vh ),
(1)
∀ vh ∈ Vh
detto Problema di Galerkin. Indicando con {ϕi ; i = 1, 2, . . . , Nh } una base di Vh ,
basta che la (1) sia soddisfatta per ogni funzione di base
Scrivendo uh =
PN h
j=1
a(uh , ϕi ) = F (ϕi ),
i = 1, 2, . . . , Nh .
uj ϕj le equazioni diventano:
Nh
X
uj a(ϕj , ϕi ), i = 1, 2, . . . , Nh .
j=1
ANALISI NUMERICA – p.2/22
METODO di GALERKIN
Si ottiene il sistema lineare
Au = f
-A ∈ RNh ×Nh : Matrice di stiffness,
Aij = a(ϕj , ϕi );
- f ∈ RNh : vettore dati, fi = F (ϕi );
- u ∈ RNh : vettore incognite, ui = ui .
Proprietà di A
A associata a problemi ellittici è definita positiva;
A simmetrica se e solo se la forma bilineare a(·, ·) è simmetrica
(Pb. con condizioni Dirichlet o miste)
Altre proprietà come condizionamento, struttura di sparsità, dipendono dalla
base considerata.
A seconda della scelta del sottospazio Vh si ha:
⇒ FEM: Metodi agli Elementi Finiti (Polinomi a tratti);
⇒ Metodi Spettrali (Polinomi algebrici definiti su tutto il dominio);
ANALISI NUMERICA – p.3/22
METODO di GALERKIN: analisi
Esistenza e Unicità
Il Lemma di Lax-Milgram vale per ogni spazio di Hilbert e quindi in
particolare per lo spazio Vh sottospazio dello spazio di Hilbert V ;
a(·, ·) e F (·) sono gli stessi;
⇒ sono soddisfatte le ipotesi del Lemma per cui si ha:
COROLLARIO
La soluzione del problema di Galerkin esiste ed è unica
Analogamente si può dimostrare l’esistenza e unicità della soluzione dal fatto
che la matrice A essendo definita positiva risulta non singolare.
Stabilità: la soluzione discreta rimane limitata al tendere di h a zero. Infatti
vale:
COROLLARIO
Il metodo di Galerkin è stabile, uniformemente rispetto ad h in quanto vale
kuh kV ≤
1
kF kV 0
α
ANALISI NUMERICA – p.4/22
METODO di GALERKIN: analisi
Si vuol dimostrare la convergenza della soluzione del problema di Galerkin alla
soluzione del problema debole quando h tende a zero.
Lemma di Céa
Il metodo di Galerkin è fortemente consistente :
a(u − uh , vh ) = 0, ∀vh ∈ Vh , Ortogonalità di Galerkin
DIM: Essendo Vh ⊂ V , la soluzione esatta u soddisfa il problema debole per ogni
v = vh :
a(u, vh ) = F (vh )
sottraendo con a(uh , vh ) = F (vh ) si ha
a(u, vh ) − a(uh , vh ) = 0
da cui la tesi per la bilinearità di a(·, ·).
ANALISI NUMERICA – p.5/22
METODO di GALERKIN: analisi
Ortogonalità di Galerkin o metodo di proiezione ortogonale:
se a(·, ·) è simmetrica allora definisce un prodotto scalare in V e la
consistenza forte dice che il vettore errore u − u h è ortogonale rispetto al
prodotto scalare a(·, ·) a tutti gli elementi di V h
La soluzione uh del metodo di Galerkin è la proiezione su V h della soluzione
esatta
⇒ minimizza la distanza dalla soluzione esatta nella norma indotta dal
prodotto scalare
u
u − uh
Vh
uh
ANALISI NUMERICA – p.6/22
METODO di GALERKIN: Convergenza
Applichiamo la forma bilineare a u − uh :
a(u − uh , u − uh )
=
a(u − uh , u − vh ) + a(u − uh , vh − uh )
=
a(u − uh , u − vh ) + 0
per la continuità di a(·, ·): |a(u − uh , u − vh )| ≤ M ku − uh kV ku − vh kV
per la coercività: a(u − uh , u − uh ) ≥ αku − uh k2V da cui mettendo insieme le
disuguaglianze trovate:
αku − uh k2V
≤
ku − uh kV
≤
a(u − uh , u − uh ) = a(u − uh , u − vh ) ≤ M ku − uh kV ku − vh kV
M
ku − vh kV , ∀ vh ∈ Vh ⇒
α
M
inf ku − wh kV
ku − uh kV ≤
α wh ∈Vh
ANALISI NUMERICA – p.7/22
METODO di GALERKIN: Convergenza
affinchè si abbia la convergenza occorre richiedere che al tendere a zero di h,
Vh tenda a V (hp di densità dei sottospazi):
lim
inf
h→0 vh ∈Vh
kv − vh kV = 0,
∀v ∈ V
ne segue che per valori limite di h (raffinando) la soluzione approssimata
converge all soluzione esatta:
lim ku − uh kV = 0
h→0
dalla disuguaglianza trovata si può studiare la velocità di convergenza (i.e.
"come" l’errore va a zero). Tale velocità dipende da come è stato scelto il
sottospazio Vh
il metodo di Galerkin genera un errore che è minore della minima distanza di v
dal sottospazio Vh moltiplicato per una costante ( M
) che non dipende da h,
α
per cui si può concludere essere ottimale.
ANALISI NUMERICA – p.8/22
ELEMENTI FINITI 1D
Ω = (a, b) ⊂ R: si vogliono delle approssimazioni di H 1 (a, b) dipendenti da un
parametro h:
sia Th una partizione di (a, b) in N + 1 sottointervalli K j = (xj−1 , xj ) di
ampiezza hj = xj − xj−1 e h = maxj hj con
a = x0 < x1 < · · · < xN +1 = b
le funzioni di H 1 (a, b) sono continue per cui possiamo costruire la seguente
famiglia di sottospazi:
Xhr = {vh ∈ C 0 (Ω) : vh |Kj ∈ Pr , ∀Kj ∈ Th }
Xhr sono sottospazi di H 1 (a, b) in quanto contengono funzioni derivabili
tranne che in un numero finito di punti (xj );
Gli elementi di Xhr sono localmente polinomi e quindi facili da calcolare;
conviene scegliere una base ϕi in modo che il supporto della generica
funzione intersechi quello di un numero piccolo di altre funzioni
(sparsità matrice rigidezza e costo computazionale ridotto)
conviene scegliere una base lagrangiana: i coefficienti di vh nella base
coincidono con i valori di vh in punti opportuni detti nodi.
ANALISI NUMERICA – p.9/22
ELEMENTI FINITI 1D: spazio
1
Xh
Spazio delle funzioni continue lineari a tratti su T h , partizione di (a, b):
gradi di libertà, ovvero i valori che occorre assegnare per individuare
univocamente un elemento dello spazio, coincidono con N + 2, il numero degli
elementi della partizione (nodi e elementi partizione coincidono)
base lagrangiana:
ϕi ∈ Xh1 t.c. ϕi (xj ) = δij

 1,
=
 0,
i=j
i 6= j
funzione lineare a tratti che vale 1 in xi e zero in tutti gli altri nodi:
ϕi (x) =





x−xi−1
,
xi −xi−1
xi+1 −x
,
xi+1 −xi



 0,
per xi−1 ≤ x ≤ xi
per xi ≤ x ≤ xi+1
altrimenti.
ANALISI NUMERICA – p.10/22
ELEMENTI FINITI 1D: spazio
1
Xh
1
0
0
0.2
0.4
0.6
0.8
1
supporto di ϕi : [xi−1 , xi+1 ] per i = 1, . . . , N ,
per i = 0, N + 1 il supporto è il solo intervallo [x 0 , x1 ] e [xN , xN +1 ]
solo le funzioni ϕi−1 , ϕi+1 hanno supporto con intersezione non nulla con
[xi−1 , xi+1 ];
⇒ la matrice di rigidezza è tridiagonale
ANALISI NUMERICA – p.11/22
ELEMENTI FINITI 1D: spazio
1
Xh
1
0
0
0.2
0.4
0.6
0.8
1
supporto di ϕi : [xi−1 , xi+1 ] per i = 1, . . . , N ,
per i = 0, N + 1 il supporto è il solo intervallo [x 0 , x1 ] e [xN , xN +1 ]
solo le funzioni ϕi−1 , ϕi+1 hanno supporto con intersezione non nulla con
[xi−1 , xi+1 ];
⇒ la matrice di rigidezza è tridiagonale
ANALISI NUMERICA – p.11/22
ELEMENTI FINITI 1D: spazio
1
Xh
le funzioni ϕi , ϕi+1 in [xi , xi+1 ], si ripetono a meno di un fattore di scalatura
legato ad hi . Nella pratica si possono ottenere trasformando due funzioni di
b0 , ϕ
b1 costruite su un intervallo di riferimento, tipicamente [0, 1]:
riferimento ϕ
(xi , xi+1 ) intervallo della partizione di (a, b) può essere ottenuto da (0, 1)
tramite:
φ : [0, 1] → [xi , xi+1 ],
le ϕ
b0 , ϕ
b1 in [0, 1]
x = φ(ξ) = xi + ξ(xi+1 − xi )
ϕ
b0 (ξ) = 1 − ξ,
le generiche ϕi , ϕi+1 in [xi , xi+1 ]:
ϕi (x) = ϕ
b0 (ξ(x)),
con ξ(x) = (x − xi )/(xi+1 − xi )
,ϕ
b1 (ξ) = ξ
ϕi+1 (x) = ϕ
b1 (ξ(x))
ANALISI NUMERICA – p.12/22
ELEMENTI FINITI 1D: spazio
2
Xh
Spazio di polinomi compositi di grado 2 su ciascun intervallo di T h ,
gradi di libertà: occorre assegnare 3 punti in ogni intervallo, per garantire la
continuità due punti coincidono con gli estremi dell’intervallo generico di T h e il
terzo sarà il punto medio, i nodi saranno: a = x0 , ..., x2N +2 = b
base lagrangiana:
ϕi ∈ Xh2 t.c. ϕi (xj ) = δij
funzione quadratica a tratti che vale 1 in x i e zero in tutti gli altri nodi:
i pari
i dispari
ϕi (x) =





(x−xi−1 )(x−xi−2 )
,
(xi −xi−1 )(xi −xi−2 )
(xi+1 −x)(xi+2 −x)
,
(xi+1 −xi )(xi+2 −xi )



 0,

 (xi+1 −x)(x−xi−1 ) ,
(xi+1 −xi )(xi −xi−1 )
ϕi (x) =
 0,
per xi−2 ≤ x ≤ xi
per xi ≤ x ≤ xi+2
altrimenti.
per xi−1 ≤ x ≤ xi+1
altrimenti.
ANALISI NUMERICA – p.13/22
ELEMENTI FINITI 1D: spazio
2
Xh
Spazio di polinomi compositi di grado 2 su ciascun intervallo di T h ,
gradi di libertà: occorre assegnare 3 punti in ogni intervallo, per garantire la
continuità due punti coincidono con gli estremi dell’intervallo generico di T h e il
terzo sarà il punto medio, i nodi saranno: a = x0 , ..., x2N +2 = b
base lagrangiana:
ϕi ∈ Xh2 t.c. ϕi (xj ) = δij
funzione quadratica a tratti che vale 1 in x i e zero in tutti gli altri nodi:
i pari
i dispari
ϕi (x) =





(x−xi−1 )(x−xi−2 )
,
(xi −xi−1 )(xi −xi−2 )
(xi+1 −x)(xi+2 −x)
,
(xi+1 −xi )(xi+2 −xi )



 0,

 (xi+1 −x)(x−xi−1 ) ,
(xi+1 −xi )(xi −xi−1 )
ϕi (x) =
 0,
per xi−2 ≤ x ≤ xi
per xi ≤ x ≤ xi+2
altrimenti.
per xi−1 ≤ x ≤ xi+1
altrimenti.
0
0
0.5
1
ANALISI NUMERICA – p.13/22
ELEMENTI FINITI 1D: spazio
2
Xh
come nel caso lineare la base può essere descritta considerando delle
funzioni di base sull’intervallo di riferimento [0, 1] e poi trasformandole nel
generico intervallo:
ϕ
b0 (ξ) = (1 − ξ)(1 − 2ξ),
ϕ
b1 (ξ) = 4(1 − ξ)ξ,
ϕ
b2 (ξ) = ξ(2ξ − 1)
0
0
0.5
per Xh2 la matrice di stiffness risulta pentadiagonale
Aumentando il grado polinomiale aumenta numero di gradi di libertà e costo
computazionale per risolvere il sistema lineare;
l’uso di gradi elevati su nodi uniformi porta ad approssimazioni meno stabili
ANALISI NUMERICA – p.14/22
In generale:
una approssimazione elementi finiti è un’ approssimazione di Galerkin in cui lo
spazio Vh :
è associato ad una triangolazione Th nel dominio Ω: suddivisione del dominio
in oggetti geometrici dello stesso tipo:
intervalli in 1D
triangoli o rettangoli in 2D
...
Su ogni elemento K di Th gli elementi di Vh sono funzioni polinomiali
Vh ammette una base ϕi , i = 1, . . . , Nh di funzioni a supporto limitato
in generale ϕi , i = 1, . . . , Nh sono polinomi a tratti di grado basso
ANALISI NUMERICA – p.15/22
Elemento Finito Lagrangiano: definizione
Un Elemento Finito si caratterizza con la terna: ( K , Σ , Pr )
K : dominio di Definizione o oggetto geometrico
Pr : spazio dei polinomi su cui è definita una base {ϕ j }rj=0 ;
i gradi di libertà aj di un generico elemento di Pr :
p ∈ Pr , p =
{αj }ri=0 :
r
X
αj ϕ j
j=0
Si indica Σ :=
la conoscenza della posizione dei nodi permette di
trovare i gradi di libertà.
b
Supponiamo di definire l’elemento a partire da un elemento di riferimento su K
e tramite una trasformazione φ si passa all’elemento finito su K . Allora se
br ) è elemento finito di Lagrange e φ : K
b , Σ
b, P
b → Rd iniettiva e continua
(K
allora ( K , Σ , Pr ) è ancora un elemento finito di Lagrange con
b
K = φ(K),
br },
Pr = {p : K → R, p ◦ φ ∈ P
b
Σ = φ(Σ)
ANALISI NUMERICA – p.16/22
Approssimazione con Elementi Finiti Lineari
PROBLEMA

 −u00 + σu = f,
 u(a) = 0,
Formulazione debole
trovare u ∈ H01 (a, b)
Z 1
Z 1
u0 v 0 dx +
σuvdx
0
a < x < b,
u(b) = 0
Z
=
0
1
0
f vdx, ∀v ∈ H01 (a, b)
Problema agli elementi finiti
sia Th una partizione di (a, b) in N + 1 sottointervalli K j , h = (b − a)/(N + 1) e
consideriamo lo spazio delle funzioni lineari a tratti nulle nel bordo
Vh = {vh ∈ Xh1 : vh (a) = vh (b) = 0} ⊂ H01 (a, b)
trovare uh ∈ Vh
Z
b
a
0
u0h vh
dx
+
Z
b
σuh vh dx
a
=
Z
b
f vh , ∀vh ∈ Vh
a
ANALISI NUMERICA – p.17/22
Approssimazione con Elementi Finiti Lineari
Consideriamo la base delle funzioni a cappuccio
imponendo le condizioni al bordo, possiamo scrivere
N
X
uh =
ui ϕ i
i=1
richiediamo che il problema sopra sia soddisfatto per ogni elemento della
base: v = ϕi
si ottiene un sistema lineare di N equazioni nelle N incognite u i = uh (xi )
A = [aij ],
aij =
u = [ui ],
Z
Au = f
b
a
ϕ0j ϕ0i dx
f = [fi ],
+
Z
b
σϕj ϕi dx
a
fi =
Z
b
f ϕi dx
a
A è tridiagonale per elementi finiti lineari;
ANALISI NUMERICA – p.18/22
Approssimazione con Elementi Finiti Lineari
Per calcolare gli elementi di A si utilizzano le funzioni definite nell’intervallo di
riferimento e poi si trasformano opportunamente gli integrali:
Z b
ϕ0j ϕ0i + σϕj ϕi dx =
aij =
Z
=
a
xi
xi−1
ϕ0j ϕ0i
+ σϕj ϕi dx +
Z
xi+1
xi
ϕ0j ϕ0i
+ σϕj ϕi dx
Consideriamo il primo addendo per j = i − 1 e applichiamo la trasformazione
x = (x(ξ)) con dx = d(xi−1 + ξhi ) = hi dξ :
Z xi
ϕ0i−1 ϕ0i + σϕi−1 ϕi dx =
xi−1
Z
1
0
ϕi−1 (x(ξ))0 ϕi (x(ξ))0 + σ(x(ξ))ϕi−1 (x(ξ))ϕi (x(ξ)) hi dξ
ϕi (x(ξ)) = ϕ
b1 (ξ) e ϕi−1 (x(ξ)) = ϕ
b0 (ξ), per le derivate si ha
dξ
1
d
1
d
ϕ
b1 (ξ)0 ,
ϕ
b0 (ξ)0
ϕi (x(ξ)) =
ϕ
b1 (ξ)0 =
ϕi−1 (x(ξ)) =
dx
dx
hi
dx
hi
Z 1
1
⇒ ai−1,i =
ϕ
b1 (ξ)0 ϕ
b0 (ξ)0 + σ ϕ
b1 (ξ)ϕ
b0 (ξ)hi dξ
h
i
0
Per coefficienti costanti gli integrali possono essere calcolati una volta per tutte.
ANALISI NUMERICA – p.19/22
Stima dell’errore in norma H 1
Errore di Interpolazione
Data v ∈ C 0 ([a, b]), chiamiamo Π1h v l’interpolante di v nello spazio Xh1 sulla
partizione Th rispetto alla base ϕi , i = 0, . . . , N + 1 lagrangiana:
Π1h v(xi ) = v(xi ),
Π1h v(x) =
N
+1
X
v(xi )ϕi (x)
i=0
TEOREMA
r
Sia v ∈ H r+1 (I), I = (a, b), r ≥ 1, Πr
h v ∈ Xh la sua interpolante. Vale la seguente stima
dell’errore di interpolazione:
|v − Πrh v|H k (I) ≤ Ck,r hr+1−k |v|H r+1 (I) ,
k = 0, 1.
- Ck,r : costante indipendente da v e h;
- H 0 (I) = L2 (I) e | · |H 0 = k · kL2
ANALISI NUMERICA – p.20/22
Stima dell’errore in norma H 1
Errore di Approssimazione del metodo degli Elementi Finiti
Grazie alla stima dell’errore di interpolazione si dimostra:
TEOREMA
Sia Ω = I , u ∈ H p+1 (I) soluzione esatta del problema variazionale e uh la sua soluzione
T
r
approssimata con il metodo degli Elementi Finiti di grado r in V h = Xh
V , r ≤ p. Vale la
seguente stima apriori dell’errore:
||u − uh ||V ≤
M
Chr |u|H r+1 (I) .
α
Per aumentare l’accuratezza si può
diminuire h, ovvero raffinare la griglia
aumentare r, ovvero utilizzare elementi finiti di grado più elevato
u ∈ V ∩ H p+1 (I), quindi il massimo valore di r che ha senso prendere è r = p:
se la soluzione non è molto regolare non conviene usare elementi finiti di
grado elevato
Se p = 0: soluzione con regolarità minima la stima non è più valida, ma si ha
comunque convergenza
ANALISI NUMERICA – p.21/22
Stima dell’errore in norma H 1
Errore di Approssimazione del metodo degli Elementi Finiti
r
p=0
p=1
p=2
p=3
p=4
1
converge
h1
h1
h1
h1
2
converge
h1
h2
h2
h2
3
converge
h1
h2
h3
h3
4
converge
h1
h2
h3
h4
In generale: se u ∈ H p+1 (I), p > 0 allora esiste una costante C indipendente da u e
h, tale che
||u − uh ||H 1 (I) ≤ Chs |u|H s+1 (I) , s = min{p, r}.
ANALISI NUMERICA – p.22/22