TECNICHE COMPUTAZIONALI AVANZATE
Francesca Pelosi e Salvatore Filippone
Università di Roma “Tor Vergata”
Problemi Ellittici
Metodo di Galerkin per Pb Ellittici: 1D, 2D
http://www.mat.uniroma2.it/∼pelosi/
TECNICHE COMPUTAZIONALI AVANZATE – p.1/59
EDP o PDE
Le Equazioni alle Derivate Parziali sono equazioni differenziali contenenti derivate
della funzione incognita rispetto a più variabili (temporali o spaziali).
Sia la funzione incognita u = u(t, x1 , . . . , xd ), una generica EDP si può esprimere
∂u ∂u
∂u
∂ p1 +···+pd +pt u
F x, t, u,
,...,
,...,
,g = 0
,
∂t ∂x1
∂xd
∂xp11 · · · ∂tpt
(1)
di ordine p = max{p1 + · · · + pd + pt }, pi ∈ N, con g insieme di dati.
Notazione:
P
operatore divergenza: ∇ · v = div(v) = di=1
operatore di Laplace (o laplaciano): ∆u =
∂vi
,
∂xi
Pd
∂2u
i=1 ∂x2
i
TECNICHE COMPUTAZIONALI AVANZATE – p.2/59
EDP o PDE
ESEMPI:
I ordine:
+ ∇ · (bu) = 0: equazione di trasporto (convezione), conservazione
della massa in un sistema materiale che occupa la regione Ω, u
rappresenta la densità del sistema e b(x, t) è la velocità posseduta da
una particella del sistema che all’istante t occupa la posizione x;
∂u
∂t
II ordine:
−∆u = f : equazione del potenziale, diffusione di un fluido in una regione
Ω omogenea ed isotropa o spostamento verticale di una membrana
elastica (ellittica)
∂u
∂t
− ∆u = 0: equazione del calore (parabolica)
∂2u
∂t2
− ∆u = 0: equazione delle onde (iperbolica)
TECNICHE COMPUTAZIONALI AVANZATE – p.3/59
Classificazione EDP del II ordine
Una EDP del II ordine a coefficienti costanti:
∂2u
∂2u
∂2u
∂u
∂u
Lu = A 2 + B
+C 2 +D
+E
+ Fu = G
∂x1
∂x1 ∂x2
∂x2
∂x1
∂x2
con G funzione assegnata e A, B, C, D, E, F ∈ R, può essere classificata a
seconda del segno del discriminante 4 = B 2 − 4AC :
se 4 < 0 : l’equazione si dice ellittica
se 4 = 0 : l’equazione si dice parabolica
se 4 > 0 : l’equazione si dice iperbolica
In generale non è possibile ricavare per via analitica una soluzione della (1):
i metodi di integrazione analitica disponibili (tecnica di separazione delle
variabili) sono di limitata applicabilità.
Anche nel caso si conosca un integrale generale non è detto si riesca a
determinare un integrale particolare.
⇒ Si ricorre a metodi numerici per ottenere una soluzione approssimata
TECNICHE COMPUTAZIONALI AVANZATE – p.4/59
Parte I
EQUAZIONI ELLITICHE
TECNICHE COMPUTAZIONALI AVANZATE – p.5/59
Eq. di Poisson
Si considera un dominio Ω ⊂ R2 limitato e connesso e sia ∂Ω la sua frontiera.
Il problema è:
−∆u = f
in Ω
dove f = f (x) è una funzione assegnata nelle due variabili spaziali x = (x 1 , x2 ).
L’equazione, ellittica del II ordine, si dice formulazione forte del problema di Poisson,
(o equazione di Laplace nel caso f = 0).
Per avere un unica soluzione occorre aggiungere opportune condizioni al contorno:
Problema di Dirichlet: si può assegnare il valore dello spostamento sul bordo:
u=g
su ∂Ω
dove g è una funzione assegnata. Se g = 0 si dice omogeneo.
Problema di Neumann: si può assegnare il valore della derivata normale di u
sul bordo:
∂u
∇u · n =
= h su ∂Ω
∂n
h è una funzione assegnata, n la normale uscente ad Ω.
Se h = 0 si dice omogeneo
TECNICHE COMPUTAZIONALI AVANZATE – p.6/59
Eq. di Poisson
Problema misto: si possono assegnare su porzioni diverse della frontiera ∂Ω
condizioni di tipo diverso.
Se ∂Ω = ΓD ∪ ΓN con ΓD ∩ ΓN = ∅ si possono imporre le condizioni:

 u = g,
 ∂u = h,
∂n
su ΓD
su ΓN .
Anche nel caso di Dirichlet omogeneo in cui f sia continua nella chiusura di Ω
(Ω) non è detto il problema ammetta soluzione regolare: u potrebbe non
appartenere a C 2 (Ω). ESEMPIO
Se f ∈
/ C 0 (Ω) non ha senso in generale cercare soluzioni u ∈ C 2 (Ω).
⇒ Siamo interessati a trovare una formulazione alternativa a quella forte che
consenta di trattare alcuni casi fisicamente significativi
TECNICHE COMPUTAZIONALI AVANZATE – p.7/59
Pb. di POISSON: 1D
Problema Dirichlet omogeneo:

 −u00 (x) = f (x),
 u(0) = 0,
0<x<1
u(1) = 0
Configurazione di equilibrio di un filo elastico con una tensione pari a 1, fissato
agli estremi e soggetto ad una forza trasversale di intensità f (x). La funzione u
descrive lo spostamento verticale del filo rispetto alla posizione a riposo u = 0.
Se si suppone di applicare un carico concentrato in uno o più punti, la
soluzione fisica esiste ed è continua ma non derivabile.
TECNICHE COMPUTAZIONALI AVANZATE – p.8/59
Pb. di POISSON: 1D
Problema Dirichlet omogeneo:

 −u00 (x) = f (x),
 u(0) = 0,
0<x<1
u(1) = 0
Configurazione di equilibrio di un filo elastico con una tensione pari a 1, fissato
agli estremi e soggetto ad una forza trasversale di intensità f (x). La funzione u
descrive lo spostamento verticale del filo rispetto alla posizione a riposo u = 0.
Se si suppone di applicare un carico concentrato in uno o più punti, la
soluzione fisica esiste ed è continua ma non derivabile.
TECNICHE COMPUTAZIONALI AVANZATE – p.8/59
Pb. di POISSON: 1D
Problema Dirichlet omogeneo:

 −u00 (x) = f (x),
 u(0) = 0,
0<x<1
u(1) = 0
Configurazione di equilibrio di un filo elastico con una tensione pari a 1, fissato
agli estremi e soggetto ad una forza trasversale di intensità f (x). La funzione u
descrive lo spostamento verticale del filo rispetto alla posizione a riposo u = 0.
Se si suppone di applicare un carico concentrato in uno o più punti, la
soluzione fisica esiste ed è continua ma non derivabile.
⇒ Serve una formulazione alternativa che consenta di ridurre l’ordine di
derivazione richiesto sulla soluzione u:
formulazione debole permette di passare da un problema differenziale del
secondo ordine ad uno del primo.
TECNICHE COMPUTAZIONALI AVANZATE – p.8/59
Pb. di POISSON (1D) : forma debole
Come si arriva ad una forma debole:
Si moltiplica l’equazione −u00 (x) = f (x), per una funzione test v (arbitraria) e
si integra in (0, 1):
00
−u v = f v ⇒ −
Z
1
00
u vdx =
0
Z
1
f vdx
0
Applicando la formula di integrazione per parti al primo integrale (per eliminare
la derivata seconda):
−
Z
1
00
u vdx =
0
Z
1
0
u0 v 0 dx − [u0 v]10
essendo u nota al bordo possiamo considerare funzioni test che si annullano
agli estremi dell’intervallo, per cui l’equazione diventa
Z
1
0 0
u v dx =
0
Z
1
f vdx.
(2)
0
TECNICHE COMPUTAZIONALI AVANZATE – p.9/59
Pb. di POISSON (1D) : forma debole
Come si definisce lo spazio delle funzioni test?
Lo spazio V delle funzioni test sarà tale che
v(0) = v(1) = 0
e poichè anche la soluzione u soddisfa tale condizione, sarà ricercata nello
spazio V .
Dobbiamo definire lo spazio in modo che tutte le operazioni introdotte abbiano
senso.
Se ad esempio u, v ∈ C 1 ([0, 1]) allora u0 , v 0 ∈ C 0 ([0, 1]) e quindi gli integrali
sono ben definiti, ma questo esclude soluzioni fisiche non derivabili
⇒ dobbiamo richiedere una regolarità minore.
TECNICHE COMPUTAZIONALI AVANZATE – p.10/59
RICHIAMI:
Alcune definizioni:
uno spazio vettoriale V con una metrica associata mediante una distanza
d:V ×V ⇒R
si dice completo se ogni successione di Cauchy risulta convergente nella distanza d;
uno spazio vettoriale V con una norma k k che risulta completo rispetto alla metrica
d(u, v) = ku − vk, ∀u, v ∈ V
si dice spazio di Banach;
uno spazio vettoriale V con una prodotto scalare <, >V che risulta completo rispetto alla metrica
d(u, v) = ku − vkV =
√
< u − v, u − v >V , ∀u, v ∈ V
si dice spazio di Hilbert;
TECNICHE COMPUTAZIONALI AVANZATE – p.11/59
RICHIAMI:
DEF: Dato uno spazio V :
si dice funzionale su V un operatore F t.c. F : V → R.
F si dice lineare se F (αu + βv) = αF (u) + βF (v), ∀u, v ∈ V , ∀α, β ∈ R
F lineare si dice limitato se ∃C > 0 t.c. |F (v)| ≤ CkvkV , ∀v ∈ V ;
si dice spazio duale di V:
V 0 = {F : V → R, F lineare e limitato} con norma
kF kV 0 =
ne segue che kF kV 0
sup
v∈V,v6=0
|F (v)|
kvkV
≤C
si dice forma su V un’ applicazione a t.c. a : V × V → R.
a si dice bilineare se è lineare rispetto ai due argomenti:
a(αu + βv, w) = αa(u, w) + βa(v, w), ∀u, v, w ∈ V , ∀α, β ∈ R
a(u, αv + βw) = αa(u, v) + βa(u, w), ∀u, v, w ∈ V , ∀α, β ∈ R
a si dice continua se ∃M > 0 t.c.
|a(u, v)| ≤ M kukV kvkV ,
a si dice simmetrica se a(u, v) = a(v, u),
∀u, v ∈ V .
∀u, v ∈ V
TECNICHE COMPUTAZIONALI AVANZATE – p.12/59
RICHIAMI: Spazi L
p
Definizione spazi a potenza p−esima integrabile, per 1 ≤ p ≤ ∞:
p
L (0, 1) =



v : (0, 1) 7→ R t.c ||v||Lp (0,1) :=
p=1:
p=2:
||v||L1 (0,1) =
||v||L2 (0,1) =
Z
Z
1
0
p
1
p
|v(x)| dx
1
0
sZ
< +∞



|v(x)|dx
1
0
|v(x)|2 dx
Proprietà: Date due funzioni ϕ, ψ : (0, 1) 7→ R;
se ϕ2 , ψ 2 sono integrabili ⇒ ϕψ è integrabile ossia
se ϕ, ψ ∈ L2 (0, 1) ⇒ ϕψ ∈ L1 (0, 1)
Si dimostra come conseguenza della disuguaglianza di Cauchy Schwarz:
Z
1
0
ϕ(x)ψ(x)dx ≤ kϕkL2 (0,1) kψkL2 (0,1)
TECNICHE COMPUTAZIONALI AVANZATE – p.13/59
RICHIAMI: Spazi di Sobolev
Spazio di Sobolev di ordine 1 in (0, 1)
H 1 (0, 1) = v ∈ L2 (0, 1) : v 0 ∈ L2 (0, 1)
Spazio di Sobolev di ordine k in (0, 1)
n
o
k
2
(j)
2
H (0, 1) = v ∈ L (0, 1) : v
∈ L (0, 1), j = 1, . . . , k
H k (0, 1), k ≥ 0, con H 0 (0, 1) = L2 (0, 1), sono spazi di Hilbert rispetto al
prodotto scalare:
< f, g >H k =
k Z
X
j=0
con norma
kf kH k
e seminorma
|f |H k =
sZ
1
f (j) g (j) dx
0
v
u k Z 1
uX
=t
(f (j) )2 dx
j=0
1
0
0
(f (k) )2 dx = kf (k) kL2 (0,1)
TECNICHE COMPUTAZIONALI AVANZATE – p.14/59
RICHIAMI: Spazi di Sobolev
Ad esempio per k = 1
< f, g >H 1 =
kf kH 1 =
sZ
Z
1
(f )2 dx +
0
|f |H 1 =
sZ
1
f gdx +
0
Z
1
(f 0 )2 dx =
0
f 0 g 0 dx
0
1
1
0
Z
q
kf k2L2 + kf 0 k2L2
(f 0 )2 dx = kf 0 kL2
Vale la Disuguaglianza di Poincaré:
kvkL2 (0,1) ≤ CΩ |v|H 1 (0,1)
∀v ∈ H 1 (0, 1)
TECNICHE COMPUTAZIONALI AVANZATE – p.15/59
RICHIAMI: Spazi di Sobolev
⇒ Affinchè
Z
1
0 0
u v dx =
0
Z
1
f vdx
0
sia ben definito la richiesta minima è che u 0 v 0 ∈ L1 (0, 1) e per la proprietà
sopra basta che u0 , v 0 siano a quadrato integrabile: u0 , v 0 ∈ L2 (0, 1).
Consideriamo lo spazio di Sobolev
1
2
0
2
H (0, 1) = v ∈ L (0, 1) : v ∈ L (0, 1)
(derivata nel senso delle distribuzioni)
Scegliamo come V il seguente sottospazio di H 1 (0, 1):
V = H01 (0, 1) = v ∈ H 1 (0, 1) : v(0) = v(1) = 0
e chiediamo f ∈ L2 (0, 1) in modo che anche
R1
0
f vdx abbia senso
TECNICHE COMPUTAZIONALI AVANZATE – p.16/59
Pb. di POISSON (1D) : forma debole
Il problema iniziale viene ricondotto al seguente problema integrale:
trovare u ∈ V :
Z
1
0 0
u v dx =
0
Z
1
f vdx
0
∀v ∈ V = H01 (0, 1).
(3)
OSS: Le funzioni di H 1 (0, 1) non sono necessariamente derivabili in senso classico,
ovvero
H 1 (0, 1) 6⊂ C 1 ([0, 1])
( es. funzioni continue a tratti stanno in H 1 (0, 1) ma non in C 1 ([0, 1])).
Quindi la soluzione può essere continua ma non derivabile.
Il problema (3) risulta equivalente ad un problema variazionale

 J(u) = min
v∈V J(v)
trovare u ∈ V :
 J(v) = 1 R 1 (v 0 )2 dx − R 1 f vdx,
2
DIM:
0
con
(4)
0
TECNICHE COMPUTAZIONALI AVANZATE – p.17/59
Pb. di POISSON (1D) : forma debole
Tornando al problema del filo elastico, se v denota lo spostamento ammissibile
dalla posizione di equilibrio u, l’equazione (3) esprime l’uguaglianza del lavoro
fatto dalle forze interne e dalle forze esterne in corrispondenza dello
spostamento v : principio dei lavori virtuali.
J(w) esprime energia potenziale globale associata alla configurazione w: un
qualsiasi spostamento dalla configurazione di equilibrio causa un incremento
dell’energia potenziale del sistema, la (4) afferma che la soluzione debole
minimizza l’energia potenziale.
TECNICHE COMPUTAZIONALI AVANZATE – p.18/59
Pb. di POISSON (1D) : forma debole
Pb. di Dirichlet non omogeneo:
u(0) = g0 ,
u(1) = g1
Ci si può ricondurre al caso omogeneo:
definiamo la funzione Rg , detta rilevamento del dato di bordo, tale che
Rg (0) = g0 ,
Rg (1) = g1 .
In questo caso possiamo prendere Rg = (1 − x)g0 + xg1 .
Se u è la soluzione del problema non omogeneo allora
u = u − Rg
è la soluzione del problema omogeneo corrispondente
TECNICHE COMPUTAZIONALI AVANZATE – p.19/59
Pb. di POISSON(1D) : forma debole

 −u00 + σu = f,
 u0 (0) = h0 ,
Pb. di Neumann:
0<x<1
u0 (1) = h1
con σ funzione positiva e h0 , h1 ∈ R.
(Se σ = 0 la soluzione non è unica, ma definita a meno di una costante additiva.)
Analogamente al Pb. di Dirichlet: si moltiplica l’equazione per una funzione
test v , si integra in (0, 1) e si applica la formula di integrazione per parti
Z 1
Z 1
Z 1
ottenendo
u0 v 0 dx +
σuvdx − [u0 v]10 =
f vdx
0
0
0
supponiamo f ∈ L2 (0, 1) e σ ∈ L∞ (0, 1), u non è nota al bordo quindi non
richiediamo funzioni test che si annullano agli estremi ma è nota la derivata:
trovare u ∈ H 1 (0, 1)
Z 1
Z 1
u0 v 0 dx +
σuvdx
0
0
=
Z
1
0
f vdx + h1 v(1) − h0 v(0),
∀v ∈ H 1 (0, 1)
OSS: Se h0 = h1 il pb. debole è caratterizzato dalla stessa equazione del caso di
Dirichlet, ma lo spazio V delle funzioni test è ora diverso.
TECNICHE COMPUTAZIONALI AVANZATE – p.20/59
Pb. di POISSON(1D): forma debole
Pb. misto omogeneo (o di Robin):

 −u00 + σu = f,
 u(0) = 0,
0<x<1
u0 (1) = 0
In questo caso si deve chiedere che le funzioni test siano nulle in x = 0
(ΓD = 0, ΓN = 1) ponendo:
HΓ1 D (0, 1) = v ∈ H 1 (0, 1) : v(0) = 0
con i soliti passaggi la formulazione debole con f ∈ L 2 (0, 1) e σ ∈ L∞ (0, 1), è
trovare u ∈ HΓ1D (0, 1)
Z 1
Z 1
u0 v 0 dx +
σuvdx
0
0
=
Z
1
f vdx,
0
∀v ∈ H 1 ΓD (0, 1)
OSS: Di nuovo si ottiene la stessa equazione del caso di Dirichlet omogeneo, ma
con diverso lo spazio V delle funzioni test.
TECNICHE COMPUTAZIONALI AVANZATE – p.21/59
Pb. di POISSON(1D): forma debole
Pb. misto (o di Robin):

 −u00 + σu = f,
 u(0) = 0,
0<x<1
u0 (1) + γu(1) = r,
γ, r > 0
Si chiede che le funzioni test siano nulle in x = 0 (Γ D = 0 ) ponendo:
HΓ1 D (0, 1) = v ∈ H 1 (0, 1) : v(0) = 0
in questo caso il termine di bordo per x = 1 derivante dall’integrazione per
parti fornisce un termine proporzionale all’incognita u:
−[u0 v]10 = −rv(1) + γu(1)v(1)
la formulazione debole risulta:
trovare u ∈ HΓ1D (0, 1)
Z 1
Z 1
u0 v 0 dx +
σuvdx + γu(1)v(1)
0
0
=
Z
1
f vdx + rv(1),
0
∀v ∈ H 1 ΓD (0, 1)
TECNICHE COMPUTAZIONALI AVANZATE – p.22/59
Pb. di POISSON (2D) omogeneo : forma debole
Consideriamo il caso di Dirichlet omogeneo, sul dominio limitato Ω ⊂ R 2 :

 −∆u = f,
 u = 0,
in Ω
su ∂Ω
Analogamente al caso monodimensionale
Si moltiplica l’equazione differenziale −∆u = f, per una funzione test v
(arbitraria) e si integra in Ω:
Z
Z
−
∆uv dΩ =
f v dΩ
Ω
Ω
Occorre applicare l’analogo multidimensionale della formula di integrazione
per parti. Si ottiene con il teorema della divergenza:
Z
Z
div(a) dΩ =
a · n dγ
Ω
∂Ω
con a(x) = (a1 (x), a2 (x))T funzione sufficientemente regolare e
n(x) = (n1 (x), n2 (x))T versore normale uscente ad Ω e x = (x1 , x2 )T vettore
coordinate spaziali.
TECNICHE COMPUTAZIONALI AVANZATE – p.23/59
Pb. di POISSON (2D) omogeneo : forma debole
Applicando il teorema della divergenza alle funzioni (0, ϕψ) T e (ϕψ, 0)T si
ottiene
Z
Z
Z
∂ϕ
∂ψ
ψ dΩ = −
ϕ
dΩ +
ϕψni dγ, i = 1, 2
∂x
∂x
i
i
Ω
Ω
∂Ω
P
∂
∂u
Sfruttando ∆u = div∇u = i=1,2 ∂x ∂x :
i
−
Z
∆uv dΩ
=
Ω
X Z
i=1,2
=
=
Ω
∂
∂xi
Z
i
∂u
∂xi
v dΩ
Z
X ∂u ∂v
X ∂u
−
dΩ −
vni dγ
Ω i=1,2 ∂xi ∂xi
∂Ω i=1,2 ∂xi
X Z ∂u ∂v
X Z ∂u dΩ −
ni v dγ
∂xi ∂xi
∂xi
i=1,2 Ω
i=1,2 ∂Ω
si arriva alla formula di Green
−
Z
∆uv dΩ =
Ω
Z
Ω
∇u · ∇v dΩ −
Z
∂Ω
∂u
v dγ
∂n
TECNICHE COMPUTAZIONALI AVANZATE – p.24/59
Pb. di POISSON (2D): forma debole
Si arriva così alla formulazione debole
trovare u ∈ H01 (Ω) :
Z
Ω
∇u · ∇v dΩ =
Z
f v dΩ,
Ω
∀v ∈ H01 (Ω).
con f ∈ L2 (Ω) e dove denotiamo
∂v
∈ L2 (Ω)}
∂xi
H 1 (Ω)
=
{v : Ω 7→ R t.c v,
H01 (Ω)
=
{v ∈ H 1 (Ω) : v = 0 su ∂Ω}
OSS: Se u, v ∈ H01 (Ω) allora ∇u, ∇v ∈ [L2 (Ω)]2 e quindi ∇u · ∇v ∈ L1 (Ω) grazie
alla disuguaglianza di Cauchy-Schwarz:
|(∇u, ∇v)| ≤ k∇ukL2 (Ω) k∇vkL2 (Ω)
TECNICHE COMPUTAZIONALI AVANZATE – p.25/59
Pb. di POISSON (2D): forma debole
Definiamo:
a : V × V → R,
F : V → R,
a(u, v) =
F (v) =
Z
Z
Ω
∇u · ∇v dΩ,
f v dΩ,
forma bilineare
funzionale lineare
Ω
Il problema debole diventa
trovare u ∈ V : a(u, v) = F (v)
∀v ∈ V.
TECNICHE COMPUTAZIONALI AVANZATE – p.26/59
Pb. di POISSON (2D) NON omogeneo
Consideriamo il seguente problema con condizioni miste non omogenee


−∆u = f,


u = g,


 ∂u
= φ,
∂n
in Ω
su ∂ΓD
su ∂ΓN
dove ΓD , ΓN costituiscono una partizione di ∂Ω, con
1
f ∈ L2 (Ω), g ∈ H 2 (ΓD ), φ ∈ L2 (ΓN ).
Grazie alla formula di Green si ottiene:
Z
e sfruttando il dato
Z
Ω
∂u
∂n
Ω
∇u · ∇v dΩ −
Z
∂Ω
∂u
v dγ =
∂n
Z
f v dΩ
Ω
= φ su ΓN :
∇u · ∇v dΩ −
Z
ΓD
∂u
v dγ −
∂n
Z
φv dγ =
ΓN
Z
f v dΩ
Ω
TECNICHE COMPUTAZIONALI AVANZATE – p.27/59
Pb. di POISSON (2D) NON omogeneo
Imponendo v = 0 in ΓD , si arriva alla seguente forma debole
trovare u ∈ Vg :
Z
∇u · ∇v dΩ
=
Ω
Z
φv dγ +
ΓN
Z
f v dΩ
Ω
∀v ∈ V
con
V
=
HΓ1 D (Ω) = {v ∈ H 1 (Ω) : v|ΓD = 0}
Vg
=
{v ∈ H 1 (Ω) : v|ΓD = g}
La formulazione sopra presenta delle problematiche:
non è simmetrica (u ∈ Vg , v ∈ V )
Vg non è un sottospazio di V .
TECNICHE COMPUTAZIONALI AVANZATE – p.28/59
Pb. di POISSON (2D) NON omogeneo
Supponiamo di conoscere una funzione continua R g , rilevamento del dato di bordo:
Rg ∈ H 1 (Ω), Rg |ΓD = g
e consideriamo
u
e = u − Rg ,
u
e|ΓD
u
e|ΓD = 0 dunque è soluzione del problema omogeneo associato e poichè
∇u = ∇e
u + ∇Rg , il problema diventa trovare u
e ∈ V = HΓ1 (Ω):
D
Z
Z
Z
Z
∇e
u · ∇v dΩ =
φv dγ +
f v dΩ −
∇Rg · ∇v dΩ ∀v ∈ V
Ω
ovvero
ΓN
Ω
Ω
trovare u
e ∈ V : a(e
u, v) = F (v)
F (v) =
Z
φv dγ
ΓN
Z
Ω
f v dΩ −
Z
Ω
∀v ∈ V
∇Rg · ∇v dΩ
Il problema è simmetrico
Le condizioni di Dirichlet sono dette essenziali in quanto imposte in modo
esplicito per determinare lo spazio di funzioni
Le condizioni di Neumann sono dette naturali in quanto soddisfatte
implicitamente dalla soluzione del problema.
TECNICHE COMPUTAZIONALI AVANZATE – p.29/59
Problemi Ellittici più Generali



Lu = f,






in Ω
u = g,
su ∂ΓD
∂u
∂nL
su ∂ΓN
= φ,
dove ΓD , ΓN costituiscono una partizione di ∂Ω, con
1
f ∈ L2 (Ω), g ∈ H 2 (ΓD ), φ ∈ L2 (ΓN ) e definendo l’operatore L:
2
X
∂
∂u
Lu = −
ai,j
+ σu
∂x
∂x
i
j
i,j=1
e la derivata conormale di u associata all’operatore L
2
X
∂u
∂u
=
ai,j
ni
∂nL
∂x
j
i,j=1
Sotto opportune ipotesi la formulazione debole rimane la stessa definendo
a(u, v)
=
F (v)
=
Z

2
X

∂u ∂v
+ σuv  dΩ
∂xj ∂xi
Ω
i,j=1
Z
Z
−a(Rg , v) +
f v dΩ +
φv dγ

ai,j
Ω
ΓN
TECNICHE COMPUTAZIONALI AVANZATE – p.30/59
Problemi Ellittici più Generali
Vediamo le ipotesi che devono essere soddisfatte:
σ(x) ∈ L∞ (Ω) e ∃α0 > 0 tale che
σ(x) ≥ α0 ,
q.o. in Ω;
i coefficienti aij : Ω → R siano funzioni continue ∀i, j = 1, 2
esista una costante α > 0 (costante di ellitticità) tale che
∀ξ = (ξ1 , ξ2 )T ∈ R2 ,
2
X
i,j=1
aij (x)ξi ξj ≥ α
2
X
ξi2 , q.o. in Ω.
i=1
TECNICHE COMPUTAZIONALI AVANZATE – p.31/59
ESEMPIO: Elasticità Lineare


−div σ(u) = fi ,


ui = 0,


 P2
j=1 σij (u)nj = gi ,
in Ω, i = 1, 2
su ΓD , i = 1, 2
su ΓN , i = 1, 2
σij (u) = λdiv(u)δij + 2µεij (u),
εij (u) =
1
2
∂uj
∂ui
+
∂xj
∂xi

=
2
X
∂σij
div σ(u) =
(u), i = 1, 2
∂xj
j=1
∂u1
∂x1
1 ∂u1
(
2 ∂x2
+
∂u2
)
∂x1
1 ∂u1
2
(
+ ∂u
)
2 ∂x2
∂x1
∂u2
∂x2
u : spostamento di un corpo elastico che occupa la regione Ω,
sotto l’azione di una forza per unità di volume di intensità f = (f 1 , f2 )T ,
soggetto ad un carico distribuito su ΓN di intensità g = (g1 , g2 )T e
fissato all’estremità ΓD .


TECNICHE COMPUTAZIONALI AVANZATE – p.32/59
ESEMPIO: Elasticità Lineare


−div σ(u) = fi ,


ui = 0,


 P2
j=1 σij (u)nj = gi ,
in Ω, i = 1, 2
su ΓD , i = 1, 2
su ΓN , i = 1, 2
con la formula di Green
Z
σ(u) : ε(v) dΩ
2 Z
X
:=
Ω
i,j=1
=
2 Z
X
i,j=1
si ottiene:
Z
Z
Γ
σ(u) : ε(v) dΩ =
Ω
λdiv(u)div(v) dΩ + 2µ
Ω
σij (u)εij (v) dΩ
Ω
Z
σij (u)nj vi dγ −
Z
Ω
f · v dΩ +
ε(u) : ε(v) dΩ =
Ω
Z
Z
2 Z
X
i,j=1
ΓN
Ω
Ω
∂σij
(u)vi dΩ
∂xj
g · v dγ
f · v dΩ +
Z
ΓN
g · v dγ
TECNICHE COMPUTAZIONALI AVANZATE – p.33/59
ESEMPIO: Elasticità Lineare in lamina con buco
una lamina sottile sottoposta a tensione uniforme lungo un asse T x = Tx ex
con un piccolo foro (r = 1) nell’interno
La distribuzione può essere alterata solo in prossimità del buco e lo stress
mantiene valori uniformi nel resto del dominio
È conveniente costruire un cerchio ipotetico con raggio uguale all’ampiezza
della lamina
Per simmetria si può risolvere il problema in un quarto di dominio
TECNICHE COMPUTAZIONALI AVANZATE – p.34/59
ESEMPIO: Elasticità Lineare


−div σ(u) = 0,


ui = p i ,


 P2
j=1 σij (u)nj = gi ,
in Ω, i = 1, 2
su ΓD , i = 1, 2
su ΓN , i = 1, 2
ΓD
ΓN
ΓN
ΓD
TECNICHE COMPUTAZIONALI AVANZATE – p.35/59
ESEMPIO: Elasticità Lineare


−div σ(u) = 0,


ui = p i ,


 P2
j=1 σij (u)nj = gi ,
Configurazione iniziale e finale;
in Ω, i = 1, 2
su ΓD , i = 1, 2
su ΓN , i = 1, 2
Distribuzione dello stress σ 1
TECNICHE COMPUTAZIONALI AVANZATE – p.36/59
ESEMPIO: Problema di Lamè


−div σ(u) = f ,


ui = 0,


 P2
j=1 σij (u)nj = gi ,
in Ω, i = 1, 2
su ΓD , i = 1, 2
su ΓN , i = 1, 2
Lamina rettangolare flessibile posta su un piano verticale, soggetta ad una forza f
perpendicolare e con un lato fissato ad un sostegno
ΓN
ΓD
ΓN
ΓN
TECNICHE COMPUTAZIONALI AVANZATE – p.37/59
ESEMPIO: Problema di Lamè


−div σ(u) = f ,


ui = 0,


 P2
j=1 σij (u)nj = gi ,
in Ω, i = 1, 2
su ΓD , i = 1, 2
su ΓN , i = 1, 2
Lamina rettangolare flessibile posta su un piano verticale, soggetta ad una forza f
perpendicolare e con un lato fissato ad un sostegno
TECNICHE COMPUTAZIONALI AVANZATE – p.37/59
Esistenza e Unicità
Lemma di Lax-Milgram
Supposto V spazio di Hilbert;
a(·, ·) : V × V → R forma bilineare coerciva:
a(v, v) ≥ αkvk2V ,
∃α > 0 t.c.
e continua:
∃M > 0 t.c.
∀v ∈ V
|a(u, v)| ≤ M kukV kvkV ,
∀u, v ∈ V
F (·) : V → R funzionale lineare limitato,
allora esiste un’unica soluzione del problema:
trovare u ∈ V : a(u, v) = F (v)
∀v ∈ V.
Corollario
La soluzione è limitata in funzione dei dati:
kukV ≤
1
kF kV 0
α
dove α è la costante di coercività associata alla forma bilineare a(·, ·), e kF kV 0 è la norma del
funzionale F.
DIM: dalla coercività di a(·, ·) si ha che αkuk2V ≤ a(u, u) = F (u) e poichè
F (u) ≤ kF kV 0 kukV , segue la tesi.
TECNICHE COMPUTAZIONALI AVANZATE – p.38/59
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 ),
∀ vh ∈ Vh
(5)
detto Problema di Galerkin. Indicando con {ϕi ; i = 1, 2, . . . , Nh } una base di Vh ,
basta che la (5) sia soddisfatta per ogni funzione di base
Scrivendo uh =
PNh
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
TECNICHE COMPUTAZIONALI AVANZATE – p.39/59
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.
TECNICHE COMPUTAZIONALI AVANZATE – p.40/59
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
α
TECNICHE COMPUTAZIONALI AVANZATE – p.41/59
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(·, ·).
TECNICHE COMPUTAZIONALI AVANZATE – p.42/59
METODO di GALERKIN: analisi
Ortogonalità di Galerkin o metodo di proiezione ortogonale:
se a(·, ·) è simmetrica allora definisce un prodotto scalare in V e la fortemente
consistenza dice che il vettore errore u − u h è ortogonale rispetto al prodotto
scalare a(·, ·) a tutti gli elemeti di Vh
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
TECNICHE COMPUTAZIONALI AVANZATE – p.43/59
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
ku − uh kV ≤
inf ku − wh kV
α wh ∈Vh
TECNICHE COMPUTAZIONALI AVANZATE – p.44/59
METODO di GALERKIN: Convergenza
affinchè si abbia la convergenza occorre richiedere che al tende a zero di h,
Vh tenda a V :
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
TECNICHE COMPUTAZIONALI AVANZATE – p.45/59
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.
TECNICHE COMPUTAZIONALI AVANZATE – p.46/59
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.
TECNICHE COMPUTAZIONALI AVANZATE – p.47/59
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
TECNICHE COMPUTAZIONALI AVANZATE – p.48/59
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
TECNICHE COMPUTAZIONALI AVANZATE – p.48/59
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))
TECNICHE COMPUTAZIONALI AVANZATE – p.49/59
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 quadratiche 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−1
per xi ≤ x ≤ xi+2
altrimenti.
per xi−1 ≤ x ≤ xi+1
altrimenti.
TECNICHE COMPUTAZIONALI AVANZATE – p.50/59
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 quadratiche 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−1
per xi ≤ x ≤ xi+2
altrimenti.
per xi−1 ≤ x ≤ xi+1
altrimenti.
0
0
0.5
1
TECNICHE COMPUTAZIONALI AVANZATE – p.50/59
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
TECNICHE COMPUTAZIONALI AVANZATE – p.51/59
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
TECNICHE COMPUTAZIONALI AVANZATE – p.52/59
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
Σ = φ(Σ)
TECNICHE COMPUTAZIONALI AVANZATE – p.53/59
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
a
f vh , ∀vh ∈ Vh
TECNICHE COMPUTAZIONALI AVANZATE – p.54/59
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;
TECNICHE COMPUTAZIONALI AVANZATE – p.55/59
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
dξ
1
d
1
ϕ
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 unaTECNICHE
volta COMPUTAZIONALI
per tutte. AVANZATE – p.56/59
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
TECNICHE COMPUTAZIONALI AVANZATE – p.57/59
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
TECNICHE COMPUTAZIONALI AVANZATE – p.58/59
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}.
TECNICHE COMPUTAZIONALI AVANZATE – p.59/59