7/23/2019 Metodi Analitici e Numerici Per l'Ingegneria Aerospaziale - Saleri
1/107
http:/
/aero.po
limi.tripod.c
om
Memory
Metodi Analitici e Numerici perlIngegneria Aerospaziale
Corso del Prof. F.E.Saleri
A.A 2006/07
7/23/2019 Metodi Analitici e Numerici Per l'Ingegneria Aerospaziale - Saleri
2/107
Questo lavoro non e stato corretto da alcun professorequasi sicuramente contiene errori
nel caso ne trovassi qualcuno scrivimi a:
questa raccolta di appunti e reperibile in forma gratuita1 al sito:
http://aero.polimi.tripod.comnessuno e autorizzato a ricavarne alcun tipo di guadagno
se non con esplicita autorizzazione dellautore.
second editionultimo aggiornamento:
8 marzo 2008
1fatto salvo per i costi di connessione
7/23/2019 Metodi Analitici e Numerici Per l'Ingegneria Aerospaziale - Saleri
3/107
Indice
1 Analisi teorica del problema modello 4
2 Metodo delle Differenze Finite 72.1 Problema modello . . . . . . . . . . . . . . . . . . . . . . . . . 7
2.1.1 Esistenza, unicita e accuratezza della soluzione . . . . . 92.2 Problema completo a coefficienti costanti . . . . . . . . . . . . 112.3 Metodo dei coefficienti indeterminati . . . . . . . . . . . . . . 122.4 Condizioni al contorno . . . . . . . . . . . . . . . . . . . . . . 14
2.4.1 Condizioni di Dirichlet . . . . . . . . . . . . . . . . . . 14
2.4.2 Risoluzione del problema di Neumann . . . . . . . . . 142.4.3 Condizioni di Robin . . . . . . . . . . . . . . . . . . . 15
2.5 Differenze Finite su intervalli non equispaziati . . . . . . . . . 162.6 Problema completo a coefficienti variabili . . . . . . . . . . . . 172.7 Caso generale . . . . . . . . . . . . . . . . . . . . . . . . . . . 182.8 Interpolazione polinomiale . . . . . . . . . . . . . . . . . . . . 192.9 Differenze Finite Compatte . . . . . . . . . . . . . . . . . . . . 20
3 Risoluzione dei sistemi lineari 233.1 Fattorizzazione . . . . . . . . . . . . . . . . . . . . . . . . . . 23
3.1.1 Gauss MEG . . . . . . . . . . . . . . . . . . . . . . . . 233.1.2 Thomas . . . . . . . . . . . . . . . . . . . . . . . . . . 24
3.2 Metodi iterativiRichardson . . . . . . . . . . . . . . . . . . . 263.2.1 Metodo del Gradiente . . . . . . . . . . . . . . . . . . 283.2.2 Metodo del Gradiente Coniugato . . . . . . . . . . . . 30
3.3 Precondizionatori . . . . . . . . . . . . . . . . . . . . . . . . . 33
4 Differenze Finite 2D e 3D 36
1
7/23/2019 Metodi Analitici e Numerici Per l'Ingegneria Aerospaziale - Saleri
4/107
INDICE
5 Spazi 42
5.1 SpazioL2 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 425.2 Spazi H1 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 425.3 Spazi X1h . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 435.4 Proprieta . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 435.5 Spazio di Hilbert . . . . . . . . . . . . . . . . . . . . . . . . . 44
6 Formulazioni deboli 456.1 Dirichlet Omogeneo . . . . . . . . . . . . . . . . . . . . . . . . 466.2 Dirichlet Non Omogeneo . . . . . . . . . . . . . . . . . . . . . 466.3 Neumann Omogeneo . . . . . . . . . . . . . . . . . . . . . . . 47
6.4 Neumann Non Omogeneo . . . . . . . . . . . . . . . . . . . . 476.5 Funzionali e Forme . . . . . . . . . . . . . . . . . . . . . . . . 486.5.1 Funzionale . . . . . . . . . . . . . . . . . . . . . . . . . 486.5.2 Forma . . . . . . . . . . . . . . . . . . . . . . . . . . . 49
6.6 Lemma di Lax-Milgram . . . . . . . . . . . . . . . . . . . . . 506.6.1 Verifica delle ipotesi . . . . . . . . . . . . . . . . . . . 51
7 Metodo degli Elementi Finiti (Galerkin) 527.1 Convergenza . . . . . . . . . . . . . . . . . . . . . . . . . . . . 55
7.1.1 Ordine di convergenza . . . . . . . . . . . . . . . . . . 557.2 Problemi di Diffusione-Trasporto(-Reazione) . . . . . . . . . . 56
7.3 Problema modello a trasporto dominante . . . . . . . . . . . . 577.3.1 Decentramento del metodo di Galerkin (UpWind) . . . 59
7.4 Problema modello a reazione dominante . . . . . . . . . . . . 607.4.1 Condensazione della matrice di massa (Mass-Lumping) 61
8 Discretizzazione di ordine 2 in spazio dipendenti dal tempo 638.1 Formulazione debole . . . . . . . . . . . . . . . . . . . . . . . 638.2 Discretizzazione in spazio con Elementi Finiti . . . . . . . . . 648.3 Discretizzazione in tempo con Differenze Finite . . . . . . . . 66
8.3.1 Eulero Esplicito . . . . . . . . . . . . . . . . . . . . . . 66
8.3.2 Eulero Implicito . . . . . . . . . . . . . . . . . . . . . . 668.3.3 Crank-Nicolson . . . . . . . . . . . . . . . . . . . . . . 668.3.4 Analisi di stabilita per il -metodo . . . . . . . . . . . 67
8.4 Discretizzazione in tempo con Elementi Finiti . . . . . . . . . 768.4.1 Elementi Finiti discontinui . . . . . . . . . . . . . . . . 87
9 Temi desame 889.1 Temi del prof. F.Saleri . . . . . . . . . . . . . . . . . . . . . . 88
9.1.1 24.Nov.05 prova in itinere I . . . . . . . . . . . . . . . 88
Appunti 2 F.Cadei BS
7/23/2019 Metodi Analitici e Numerici Per l'Ingegneria Aerospaziale - Saleri
5/107
INDICE
9.1.2 03.Feb.06 prova in itinere II . . . . . . . . . . . . . . . 899.1.3 22.Feb.06 appello desame . . . . . . . . . . . . . . . . 919.1.4 05.Lug.06 appello desame . . . . . . . . . . . . . . . . 929.1.5 01.Set.06 appello desame . . . . . . . . . . . . . . . . 939.1.6 08.Feb.07 prova in itinere II . . . . . . . . . . . . . . . 95
9.2 Tema del prof. P.Zunino . . . . . . . . . . . . . . . . . . . . . 979.2.1 19.Feb.08 appello desame . . . . . . . . . . . . . . . . 97
Appunti 3 F.Cadei BS
7/23/2019 Metodi Analitici e Numerici Per l'Ingegneria Aerospaziale - Saleri
6/107
http:/
/aero.po
limi.tripod.c
om
Capitolo1Analisi teorica del problema modello
u =f in(0, 1)u(0) =u(1) = 0
questo problema e adimensionale e omogeneo, con il problema omogeneoho una retta come soluzioneu = 0, u(0) = u0, u(1) = u1. Integrandolequazione si ottiene
u(x) = c0+c
1x
x
0
F(s) ds
con F(s) =
s0
f(t) dt s0
F(s) ds = integro per parti
=
sF(s)x
0 x
0
sF(s) ds
= xF(x) 0 x
0
sf(s) ds
= x x0
f(s) ds x0
sf(s) ds s0
F(s) ds =
x0
(x s)f(s) ds
calcolo le costanti c0 e c1
u(0) = 0 c0 = 0 perche 0
0
4
7/23/2019 Metodi Analitici e Numerici Per l'Ingegneria Aerospaziale - Saleri
7/107
1. Analisi teorica del problema modello
u(1) = 0 c1 10 (1 s)f(s) ds= 0c1 =
10
(1 s)f(s) ds
u(x) = x
10
(1 s)f(s) ds x
0
(x s)f(s) ds
=
x0
(x xs x+s)f(s) ds+ 1
x
x(1 s)f(s) ds
= x0
s(1 x)f(s) ds+ 1x
x(1 s)f(s) ds
introduco la funzione di Green
G(x.z) =
s(1 x) se s (0, x)x(1 s) se s (x, 1)
u(x) =1
0G(x, s)f(s) ds
Figura 1.1: Funzione di Green.
G e simmetrica, positiva, continua posso integrarla sempre, quindiu(x)
se f e integrabile f continua C0
se f C0(0, 1) u C2
mau =f se f e continua u C2, se f Cm u Cm+2, sempre piuregolare.Stima a priori, cerco la dipendenza continua dai dati, ossia se perturbo non
Appunti 5 F.Cadei BS
7/23/2019 Metodi Analitici e Numerici Per l'Ingegneria Aerospaziale - Saleri
8/107
7/23/2019 Metodi Analitici e Numerici Per l'Ingegneria Aerospaziale - Saleri
9/107
http:/
/aero.po
limi.tripod.c
om
Capitolo2Metodo delle Differenze Finite
Nei metodi numerici per la risoluzione diu =finvece di calcolare linte-grale si utilizza una sommatoria sui nodi di quadratura, sommatoria pesata.
k
wkG(xk, sk)f(sk) 1
0
G(x, s)f(s) ds
questo approccio viene trascurato perche lungo; si utilizzano le DifferenzeFinite.
2.1 Problema modello
u =f x (0, 1)la derivata e un operazione di passaggio al limite:
- 1passo: discretizzo, non ho un continuo da 0 a 1 ma un insieme dipunti
x0= 0 xN+1= 1
u(xi) =f(xi) i= 1 . . . N vengono esclusi gli estremi perche sono fissate le condizioni al contorno,semplificando xk =kh(distribuzione di nodi uniforme) h=
1N
- 2passo: approssimo,
D2ui = f(xi)con D2ui u(xi)
u(xi) =
u(xi)u(xi1)
h approx decentrata indietro (o avanti)
u(xi+1)u(xi1)2h
approx centrata
7
7/23/2019 Metodi Analitici e Numerici Per l'Ingegneria Aerospaziale - Saleri
10/107
2. Metodo delle Differenze Finite
si preferisce lapprossimazione centrata perche piu accurata
u(xi1) = u(xi) hu(xi) + h2
2u(xi) +. . .
u(xi+1) = u(xi) +hu(xi) +
h2
2u(xi) +. . .
sommo e risolvo rispetto au(xi)
u(xi) = 1
h2
u(xi+1) 2u(xi) +u(xi1)
+oh2 . . .
sostituendo lapprossimazione il problema diventa
1h2ui+1 2ui+ui1= f(xi)u0=uN+1= 0cerco una funzione{uk}N+1k=0 t.c. ui=u(xi)1 con ui h0 u(xi)
u0= 0
1h2
(u0 2u1+u2) =f(x1) 1
h2(u1 2u2+u3) =f(x2)
. . .
1h2
(uN1 2uN+uN+1) =f(xN)uN+1= 0
si ottiene un sistema lineare
u =
uiN+1
i=0
f =
0, f(x1), . . . , f (xN), 0
A= 1
h
2
h2 0 01 2 1 ...
... 1 2 1 ...
.
.. .
. . .
. . .
. .
.
..... 1 2 10 0 h2
si ottiene un metodo approssimato, Au= f, si puo anche ridurre alle incog-nite interne preoccupandosi di trattare opportunamente i termini noti. Conquesto metodo non si ottiene un funzione vera e propria, ma tanti punti, u euna soluzione nodale.
1ui e unapprossimazione di u(xi),u(xi) e il valore puntuale della soluzione, uiil valorenello stesso punto dellapprossimazione
Appunti 8 F.Cadei BS
7/23/2019 Metodi Analitici e Numerici Per l'Ingegneria Aerospaziale - Saleri
11/107
2. Metodo delle Differenze Finite
2.1.1 Esistenza, unicita e accuratezza della soluzione
!uh se Anon e singolare.A e simmetrica e definita positiva2
xTAx 0 x RN= 0 x= 0
verificoxTAx 0
xTAx = [x1, x2, . . . , xN1, xN]
2 1 0
1 2
1
...... . . . . . . . . . ...... 1 2 10 1 2
x1x2...
xN1xN
= [x1, x2, . . . , xN1, xN]
2x1 x2x1+ 2x2 x3
...xN2+ 2xN1 xN
xN1+ 2xN
= 2x21 x2x2 x1x2+ 2x22 x2x3 . . . xN1xN+ 2x2N
= x21+x21 2x1x2+x22+x22 . . . 2xN1xN+x2N+x2N
= x21+ (x1 x2)2 + (x2 x3)2 +. . .+ (xN1 xN)2 +x2Ne sicuramente 0, qualsiasi vettore prendo e puo essere = 0 solo nel casoin cui il vettore sia identicamente nullo, grazie al primo e allultimo termine A e definita positiva; A e invertibile, A e non singolare, !uh, esiste lapprossimazione alle Differenze Finite (ottimo, il problema inizialeaveva una soluzione e la soluzione approssimata anche).Dallapprossimazione nodale: uh(xi) = ui dove uh(xi) e la funzione dis-
cretizzata calcolata in xi mentre ui e la funzione approssimata; bisogna ac-contentarsi di confrontare lapprossimazione con la funzione calcolata neinodi.
maxi=1,...,N
|u(xi) uh(xi)|= max
i=1,...,N|eh(xi)| ehh,
ehh, h0 02autovalori con parte reale positiva
Appunti 9 F.Cadei BS
7/23/2019 Metodi Analitici e Numerici Per l'Ingegneria Aerospaziale - Saleri
12/107
2. Metodo delle Differenze Finite
il metodo e convergente, lerrore tende a zero.
metodo e convergenteconsistente e stabilecalcolo lerrore di consistenza3, il metodo e consistente se lerrore di consis-tenza tende a zero quando il parametro di discretizzazione (h) va a zero
hh0 0
h(xi) =
u(xi+1) 2u(xi) +u(xi1)
h2
f(xi)
h= maxi=1,...,N
|h(xi)| limitatoresviluppo in serie in xi
= 1h2
u(xi) +hu
(xi) +h2
2u(xi) +
h3
6u(xi) +
h4
24uv() +
2u(xi) ++u(xi) hu(xi) + h
2
2u(xi) h
3
6u(xi) +
h4
24uv()
f(xi)
(xi; xi+1)
(xi1; xi)
semplificando
= u(xi) h2
24
uv() +uv()
f(xi)
u(xi) f(xi) = 0 perche e il problema iniziale= h
2
24
uv() +uv()
xi1 < < xi+1 per il th. della media
= h2
242uv(i)
h(xi) = h2
12uv(i)
uv(i) non lo conosco e ne ho uno per intervallo
u =fuv =f
= h2
12f(i)
3e lerrore che si commette chiedendo alla soluzione esatta di soddisfare lo schemanumerico
Appunti 10 F.Cadei BS
7/23/2019 Metodi Analitici e Numerici Per l'Ingegneria Aerospaziale - Saleri
13/107
2. Metodo delle Differenze Finite
h(xi) =
h2
12 maxi=1,...,N |f
(i)|meglio limitare con il massimo di tutti
h = h2
12maxf
hh0 0 con ordine 2
se: u C4 oppure f C2uhh, 18f
la soluzione discreta e controllata dai dati4.
Th. di convergenza
Se f C2 uh h0 unei nodi di discretizzazione e
u uh h2
96
812
f
2.2 Problema completo a coefficienti costantiu +u +u = f in(0, 1)u(0) =u(1) = 0
sostituisco (0,1) con una discretizzazione di passo h (vera solo nei nodi)
u (xi) +u(xi) +u(xi) =f(xi) i= 1, . . . , N
ora si sostituiscono le derivate con i rapporti incrementali, utilizzando rappor-ti incrementali centrati anche per la derivata prima mantengo unaccuratezzadi convergenza sullerrore del secondordine
ui+1 2ui+ui1h2
+ui+1 ui1
2h +ui= f(xi) i= 1, . . . , N
4simile a: Ax= bx= A1bx A1b
Appunti 11 F.Cadei BS
7/23/2019 Metodi Analitici e Numerici Per l'Ingegneria Aerospaziale - Saleri
14/107
2. Metodo delle Differenze Finite
la matrice Arisulta essere0
. . . . . . . . .
h2
2h
2h2
+ h2
+ 2h
. . . . . .
. . .
0
non e simmetrica.
2.3 Metodo dei coefficienti indeterminati
Il metodo dei coefficienti indeterminati e un modo per costruire rapportiincrementali per differenze finite utilizzando sviluppi in serie relativi al nododinteresse e con la possibilita di fissare lordine di convergenza dellerrore.
u(xi) a1
u(xi1) +a0
u(xi) +a1
u(xi+1)
12h
0 1
2hoppure
1
h
1
h 0
ne segue
u(xi) q
j=p
aiju(xij )
si puo notare che questo metodo non dipende solo dai punti vicini al puntodi interesse, ma e estendibile a punti diversi e a piu punti.Il numero di nodi e la loro posizione sono gradi di libert a che vengono sceltiin modo da massimizzaare lordine dellapprossimazione; per es. si ricavalapprossimazione di ordine massimo della derivata prima in xi con p= 1 e
q= 1 (3 nodi)u(xi) =a1u(xi1) +a0u(xi) +a1u(xi+1)
lerrore di discretizzazione in questo caso e errore di troncamento
h(xi) = u(xi)
a1
u(xi) hu(xi) + h
2
2u(xi) h
3
6u(xi) +. . .
+a0u(xi)
+a1
u(xi) +hu
(xi) +h2
2u(xi) +
h3
6u(xi) +. . .
Appunti 12 F.Cadei BS
7/23/2019 Metodi Analitici e Numerici Per l'Ingegneria Aerospaziale - Saleri
15/107
2. Metodo delle Differenze Finite
raccogliendo i termini simili
h(xi) = u(xi)( a1 a0 a1 =0 condizione minimale
) +
+u(xi)(1 +ha1 ha1 =1 2a condizione
) +
+u(xi)(h2
2a1 h
2
2a1
=0 3a condizione
) +
+u(xi)(h3
6
a1
h3
6
a1 termine residuo
) +. . .
si intende minimizzare lerrore relativo allapprossimazione della derivata pri-ma, le condizioni da imporre sono quindi: annullare il termine inu(xi), porrela 2a condizione= 0 e annullare lerrore dovuto al termine di derivata prima,la 3a condizione = 0 serve ad aumentare lordine di convergenza dellerrore.Non si possono imporre altre condizioni perche il supporto e basato su 3 nodi,si intuisce che per aumentare lordine di convergenza e necessario allargare ilsupporto a piu nodi per poter imporre piu condizioni, ma questo da problemiai bordi del dominio. Il sistema ottenuto e
a1+a0+a1 = 0
ha1 ha1= 1h2
2a1 h22a1 = 0
a1= a12ha1= 1
a1=
12h
a0= 0
a1 = 12hrisulta uno schema centrato
u(xi) =ui+1 ui1
2h
lerrore di troncamento risulta
h(xi) u(xi)1h
h3
6h(xi) oh2
lerrore converge con ordine 2.
Appunti 13 F.Cadei BS
7/23/2019 Metodi Analitici e Numerici Per l'Ingegneria Aerospaziale - Saleri
16/107
2. Metodo delle Differenze Finite
2.4 Condizioni al contorno
2.4.1 Condizioni di Dirichlet
Le condizioni al contorno di Dirichlet fissano i valori della funzione incognitanegli estremi
u =fu(0) =u0
u(1) =u1
come gia visto una volta scritta la matrice del sistema la prima e lultima
equazione possono passare al termine noto riducendo il problema ai nodiinterni.
2.4.2 Risoluzione del problema di Neumann
Le condizioni di Neumann fissano i valori delle derivate della funzione uneinodi estremi, se applicate al problema modello generano infinite soluzioni inquanto si ottiene una soluzione a meno di una costante che a questo puntopuo assumere qualsiasi valore. Solitamente queste condizioni vengono appli-cate ad esempio al problemau +u =foppure miste come nel problemaDirichlet-Neumann. In questo caso si e interessati al metodo di risoluzione
piu che alla soluzione stessa, verra utilizzato quindi il problema modello concondizioni al contorno miste del tipo Dirichlet-Neumann, la procedura restacomunque completamente estendibile agli altri casi
u =fu(0) =u0
u(1) =1
discretizzando si ottiene
ui+12ui+ui1
h2 =f(xi) i= 1, . . . , N u(0) =u0
u(1) approx opportuna= 1si utilizza il metodo dei coefficienti indeterminati per approssimare la derivataprima ed ottenere uno schema accuarato almeno di ordine 2; e necessariocalcolare uno schema decentrato in quanto quello centrato e si accurato diordine 2 ma non applicabile agli estremi poiche si richiederebbe lesistenza
Appunti 14 F.Cadei BS
7/23/2019 Metodi Analitici e Numerici Per l'Ingegneria Aerospaziale - Saleri
17/107
2. Metodo delle Differenze Finite
del nodo xN+2 che non esiste
u(xN+1) a0uN+1+a1uN+a2uN1 DuN+1ai2
i=0t.c. DuN+1 u(xN+1) h0
h20
sviluppo in serie
DuN+1 = a0uN+1+
+a1
uN+1 huN+1+
h2
2uN+1
h3
6u()
+ (xN, xN+1)
+a2uN+1 2huN+1+ 4 h22 uN+1 8h36 u()+ (xN1, xN)u(xN+1) DuN+1 = uN+1(a0+a1+a2) +
+uN+1(1 +a1h+ 2a2h) +
uN+1h2
2(a1+ 4a2) + o(h
3ak) termine dellerrore
il termine relativo allerrore di troncamento mostra che al massimo lerrore
puo convergere con ordine 3 rispetto ad h, non va trascurato che esiste untermine ak legato ai valori dei coefficienti da determinare che puo abbassaretale ordine
a0+a1+a2 = 0
a1h+ 2a2h= 1a1+ 4a2= 0
a2=
12h
a1= 4a2
a0 =
32h
a2 = 12h
a1 = 2hu e unapprossimazione di ordine 2: o( h
3
h) o(h2)
u(xN+1
)
3
2hu
N+1 2
hu
N+ 1
2hu
N1ui+12ui+ui1
2h =f
u(0) =u03
2huN+1 2h uN+ 12h uN1=1
2.4.3 Condizioni di Robin
Le condizioni al contorno di Robin si possono vedere come il caso generaleda cui si possono ricavare le condizioni di Dirichlet e Neumann, sono della
Appunti 15 F.Cadei BS
7/23/2019 Metodi Analitici e Numerici Per l'Ingegneria Aerospaziale - Saleri
18/107
2. Metodo delle Differenze Finite
forma u =f su(0, 1)0u(0) +0u
(0) =0
1u(1) +1u(1) =1
con 0, 0, 1, 1 R assegnati, considerando
i= 0 Dirichlet
i= 0 Neumann
2.5 Differenze Finite su intervalli non equis-paziati
Lidea e di realizzare uno schema di calcolo che risolva piu dettagliatamentecerte zone del dominio utilizzando pochi nodi dove non si e interessati amaggior dettaglio. Il problema classico e la risoluzione dello strato limite
Figura 2.1: Strato limite e ipotesi di discretizzazione.
Figura 2.2: Esempio nodi non equispaziati.
si intende ricavare uno schema che approssimi la derivata prima in xi, possi-
Appunti 16 F.Cadei BS
7/23/2019 Metodi Analitici e Numerici Per l'Ingegneria Aerospaziale - Saleri
19/107
2. Metodo delle Differenze Finite
bilmente del secondo ordine
u(xi+1) = u(xi) +hiu(xi) +
h2i2
u(xi) +h3i
6u(xi) +. . .
u(xi1) = u(xi) hi1u(xi) +h2i1
2 u(xi)
h3i16
u(xi) +. . .
sottraendo la prima equazione dalla seconda
u(xi+1) u(xi1) = (hi+hi1)u(xi) +12
u(xi)(h2i h2i1) +. . .
u(xi) u(xi+1) ui1hi+hi1
il metodo e accurato di ordine 1 rispetto adh, infatti lerrore di troncamentorisulta
h =1
2()(hi hi1
h2ih2
i1hihi1
)
h o(h) con h= max(h).Anche la derivata seconda va di ordine 1, per ottenere uno schema accuratodi ordine 2 e necessario passare ad un supporto a 5 nodi con uno schemacentrato
supp. UNIFORME supp. NON UNIFORME3 nodio(h2) 5 nodio(h2)
2.6 Problema completo a coefficienti variabili
5 (x)u diffusione
+
(x)u
trasporto
+(x)u reazione
=f(x)
si consideri il problema ridotto piu semplice
(x)u =f(x)si operi una sostituzione
w(x) =(x)u(x)
w(x) =f(x)F.D.
centr. unif.
wi=(xi)
ui+1ui12h
wi+1wi12h
=f(xi)
5e necessario prestare attenzione al senso fisico dellerquazione;u u potrebbenon risolvere il problema dinteresse
Appunti 17 F.Cadei BS
7/23/2019 Metodi Analitici e Numerici Per l'Ingegneria Aerospaziale - Saleri
20/107
2. Metodo delle Differenze Finite
lo schema e applicabile a tutti i nodi interni, quindi per i= 1, . . . , N risostituendo
14h2
(xi+1)(ui+2 ui) (xi1)(ui ui2)
= f(xi)
si giunge a un metodo a 5 punti non applicabile ai nodi i = 1e N1, inoltrese diventa costante non si ottiene piu il vecchio schema, ce qual quadrache non cosa.Osservando lo schema ottenuto si puo notare che la funzione u compare solonei nodi i 2, i , i+ 2, si puo quindi pensare di ridefinire wi su nodi fittizi
wi= (xi)ui+12 ui
12
h
la funzione u non esiste se non nei nodi, sostituendo, pero, lo schema che siottiene e
w(x) = f(x)w(xi) =
wi+ 12 wi 12h
(ui+ 12
ui) (ui ui 12 )h2
= f(xi)
(xi+ 12
)(ui+1 ui) (xi 12 )(ui ui1)h2
= f(xi)
la funzione e valutata in xi 12
, ma e una funzione nota e valutabileovunque, la funzione u e valutata su 3 nodi, lordine di accuratezza ottenutoe o(h2) e per = cost si ottiene il vecchio schema (e stato presentato untrucco per controllare lampiezza del supporto).
2.7 Caso generale
u +u flusso
+u = f
come primaw= u +uw +u = f
wi = (xi)ui+12
ui 12
h +(xi)ui
wi+12
wi 12
h +(xi)ui=f(xi)
Appunti 18 F.Cadei BS
7/23/2019 Metodi Analitici e Numerici Per l'Ingegneria Aerospaziale - Saleri
21/107
2. Metodo delle Differenze Finite
si consideri il termine di trasporto, in quanto il termine di diffusione e giastato analizzato precedentemente
+(xi+ 1
2)ui+ 1
2 (xi 1
2)ui 1
2
h + =f(x)
ue una funzione nodale, ne segue che non esiste sui nodi fittizi, il problema sirisolve approssimando il valore nel nodo fittizio con la media dei valori dellafunzione nei nodi adiacenti
ui+ 12 ui+ui+1
2
viene utilizzata la media come valore approssimato della funzione nel nodo
fittizio, sediventa costante si ritorna allo schema vecchio senza perturbarelordine del sistema; va prestata attenzione al tipo di media da impiegare,in quanto il problema potrebbe avere una direzione preferenziale, nel caso siutilizzeranno dei pesi per il calcolo della media.La stringa completa di un nodo interno risulta:
. . .i 12
2h
i12h2
;i+12
+i 12h2
+i 12
i+ 122h
+i; i 12
2h
i 12h2
. . .
2.8 Interpolazione polinomialeSi creino xi nodi con i = 0, 1, . . . , n, f(xi) e la valutazione della funzionecontinua nei nodi di discretizzazione, si costruisca un polinomio interpolatore
n f Pn tale per cui il suo valore nei nodi sia esattof(xi) =
n
f(xi)
lerrore di interpolazione risulta
err = f(x) n
f(x) =
fn+1()
(n+ 1)!
n
k=0
(x xk)se f Cn+1, massimizzo
maxx
|f(x)
n
f(x)| cnhn+1 maxx
|f(x)n+1|
Df(xi) =
n
f
(xi) i= 0, . . . , n
D2f(xi) =
n
f
(xi)
Appunti 19 F.Cadei BS
7/23/2019 Metodi Analitici e Numerici Per l'Ingegneria Aerospaziale - Saleri
22/107
2. Metodo delle Differenze Finite
scritto sulla base di Lagrange il polinomio interpolatore esiste unico
n
f(x) =n
i=0
f(xi)i(x)
Pn t.c. i(xi) =y =
1 se i= j
0 altrimenti
i =n
k=0=j
x xkxj xk polinomi caratteristici di Lagrange
Df(x) =
n
i=0 f(xi)i(x)se valuto nel punto
Df(xk) =n
i=0
f(xi)i(xk)
da questo approccio si ottengono due pregi
1. si ottengono i coefficienti ovunque senza preoccuparsi del bordo
2. maxx |f(x) n f(x)| cnhnmax|fn+1(x)| lerrore e noto (anche se
decade di un ordine rispetto ad h in quanto derivando si divide per h)i difetti
1. si ottengono approssimazioni globali, il supporto e generato dal gradoche sto usando, non si sceglie il supporto
2. su griglie equispaziate allaumentare del grado non migliora lapprossi-mazione (oscillazioni tra i nodi)
rimangono comunque degli ottimi metodi spettrali, ma non usando griglieequispaziate.
2.9 Differenze Finite Compatte
Solitamente si fissa il supporto e si pretende che
Dfi =1
j=1
ajf(xi+j )
non si considerano i valori della funzione nei punti vicini.
Appunti 20 F.Cadei BS
7/23/2019 Metodi Analitici e Numerici Per l'Ingegneria Aerospaziale - Saleri
23/107
2. Metodo delle Differenze Finite
Si consideri il problemau =f
e lo si discretizzi nel seguente modoi
Dui=
i
Sfi =
i
Sui=
i
Dui
In questo metodo si cerca di ricavare i coefficienti (pesi) direttamente mini-mizzando lerrore di troncamento
h(xi) =1
j=1 bj f(xi+j) sol. vera in schema num.
1
j=1 ajf(xi+j)= b1 f
(xi1) f(xi)hf(xi)+
h2
2 f(xi)+...
+b0f(xi) +b1 f
(xi+1) f(xi)+hf(xi)+
h2
2 f(xi)+...
= a1 f(xi1) f(xi)hf(xi)+
h2
2 f(xi)+...
a0f(xi) a1 f(xi+1) f(xi)+hf(xi)+
h2
2 f(xi)+...
(2.1)
il numero di termini da utilizzare nello sviluppo in serie dipende dal sup-porto che si intende utilizzare, nel caso di supporto a tre nodi si generano6 coefficienti, si possono scrivere fino a 6 equazioni; in realt a si utilizzano 5equazioni tutte in funzione di uno dei coefficienti in quanto con 6 equazionisi ottiene solo la soluzione banale. In questo caso quindi la serie va interrottaal termine di derivata quarta, le equazioni si scrivono intese ad annullarelerrore di troncamento
f(xi) a1+a1+a0 = 0
f(xi) b1+b1+b0+ha1 ha1= 0f(xi) hb1+hb1 h22a1 h
2
2a1= 0
f
(xi)
h2
2b1+
h2
2b1+
h3
6a1 h3
6a1= 0fv(xi) h36b1+ h3
6b1 h424 a1 h
4
24a1= 0
a1+a1+a0= 0
b1+b1+b0+ha1 ha1= 02b1+ 2b1 ha1 ha1= 03b1+ 3b1+ha1 ha1= 04b1+ 4b1 ha1 ha1= 0
b1 = b1
b0 = 4b1
a1 = 3b1
h
a1= 3b1ha0 = 0
Appunti 21 F.Cadei BS
7/23/2019 Metodi Analitici e Numerici Per l'Ingegneria Aerospaziale - Saleri
24/107
2. Metodo delle Differenze Finite
b1Dfi1+b0Dfi+b1Dfi+1 = a1f(xi1) +a0f(xi) +a1f(xi+1)
b1Dfi1+ 4b1Dfi+b1Dfi+1 = 3b1h
f(xi1) +3b1
h f(xi+1)
se b1= 0 semplificoDfi1+ 4Dfi+Dfi+1 =
3
h
f(xi+1) f(xi1)
si ottiene uno schema su 3 nodi come desiderato, lordine di convergenzadellerrore risulta essere
h fvh
4
24b1+
h4
24b1 =0, b1=b1
+h5
120a1 h
5
120a1
fv
3b1h
h5
1203b1
h
h5
120
b1
40fv
h5
hh o(h4)
si deduce che si ottiene accuratezza di ordine 4 rispetto ad h con un supportodi 3 nodi.Il difetto di questo metodo e il non poter conoscere il valore della derivata in
un punto, da informazioni globali; inoltre come si generano Dfi1, Dfi, Dfi+1?Per il problema modello
u = fu Suf M u
Su = MuSf = Mu
risulta essere tutto noto eccezion fatta per la funzione u incognita.
Appunti 22 F.Cadei BS
7/23/2019 Metodi Analitici e Numerici Per l'Ingegneria Aerospaziale - Saleri
25/107
http:/
/aero.po
limi.tripod.c
om
Capitolo3Risoluzione dei sistemi lineari
Con i metodi appena visti si giunge a formulare il problema nella forma
Ax= b
e necessario quindi ottimizzare il processo di risoluzione, A e una matricetri-diagonale
A=
a11 a12 0
a21 a22. . .
.. .
.. .
.. .. . . . . . an1n
0 ann1 ann
per risolvere il sistema in x e necessario invertire la matrice A, risulta piucomodo fattorizzarla e risolvere due sistemi tridiagonali.
3.1 Fattorizzazione
3.1.1 Gauss MEG
La fattorizzazione di Gauss1 consiste nello scomporre la matrice di partenzain due matrici triangolari tali per cui sia verificato
A= LU=
1 0 1...
. . .
. . . 1
. . .
...
. . . 0
1Calcolo Numerico
23
7/23/2019 Metodi Analitici e Numerici Per l'Ingegneria Aerospaziale - Saleri
26/107
3. Risoluzione dei sistemi lineari
il metodo di eliminazione di Gauss (MEG) perviene alla matrice triangolareeliminando una colonna alla volta
Figura 3.1: Schema di fattorizzazione.
a(k+1)ij =a
(k)ij mija(k)ij i, j =k + 1 . . . n
mik = a(k)ik
a(k)kk pivot
se akk= 0lidea e quella di, anziche invertire la matriceAper risolvere il sistema lineare,risolvere due sistemi triangolari, computazionalmente meno oneroso
Ly=bUx= y
3.1.2 Thomas
Nel caso particolare di A tri-diagonale si puo utilizzare una variante delmetodo di Gauss, lalgoritmo di Thomas, la fattorizzazione produrra matricidel tipo
L=
1 0
l21 1. . . . . .
0 lnn1 1
, U= u11 u12 0
u22 . . .. . . un1n
0 unn
facendo il prodotto
1ar 1ac u11 = a112ar 2ac u12 = a122ar 4ac u23 = a23
Appunti 24 F.Cadei BS
7/23/2019 Metodi Analitici e Numerici Per l'Ingegneria Aerospaziale - Saleri
27/107
3. Risoluzione dei sistemi lineari
tutti gli elementi della sopradiagonale diUsono uguali a quella di A
2ar 1ac l21 u11=a21l21=
a21u11
=a21a11 noto da prima
2ar 2ac l21 a12+u22 =a22u22 = a22 a21a12
a113ar 3ac l32 u22=a32
l32=a32u22
li+1,i =
ai+1,i
uiinoto a priori
noto dal passo precendente
il costo computazionale e o(n).Va ricordato che nel caso di matrici simmetriche e definite positive e parti-colarmente efficiente la fattorizzazione di Cholesky
A= RTR con R triangolare superiore
In realta si sta risolvendo un sistema perturbato (60 Wilkinson)
P,L,U= lu(A)con Psi intende la matrice delle permutazioni
A+A
x+x x
= b+b
si cerca di capire da cosa dipende la propagazione degli errori
A = 0 x x
x K(A)bb
con K(A) numero di condizionamento2 della matrice A,
K(A) = AA1se A e s.d.p.
K(A) =max(A)
min(A)
2e un indice di quanto gli errori di approssimazione introdotti per effetto delle ap-prossimazioni durante i calcoli influiscano sulla validita della soluzione trovata risolvendoil sistema lineare
Appunti 25 F.Cadei BS
7/23/2019 Metodi Analitici e Numerici Per l'Ingegneria Aerospaziale - Saleri
28/107
3. Risoluzione dei sistemi lineari
la matrice risulta ben condizionata se K(A) e vicino a 1.
A2= supx=0
Ax2x2 = maxx=1 Ax2
conATAraggio spettrale (autovalore massimo); nei casi di interesse, matricigenerate da metodi alle Differenze Finite si ottiene
A= 1
h2
2 1 0
1 . . . . . .. . . . . .
1
0 1 2
K(A) 1
h2 K(A) h0
riducendo il lintervallo h tra i nodi miglioro la discretizzazione e riducolerrore di troncamento ma si uccide la veridicita della soluzione trovata. Perrisolvere questo problema o si complica il metodo di risoluzione o si cambiail metodo stesso.
3.2 Metodi iterativi Richardson
Ax = bA = P N splitting additivo
P x = b+N x
dato x(0) si calcola per k= 0, 1, . . .
x(k+1) =P1(b+N x(k))
ci si aspetta x(0) arbitrario, il metodo e computazionalmente poco oneroso,convergente lim
kx(k) = x. Pva scelta facile da invertire per alleggerire il
metodo, vine scelta diagonale, tri-diagonale, triangolare (e lunica matriceche va invertita), ottimo se P =I
e(k) = x(k) xP(x x(k+1)) = N(x x(k))
e(k+1) = P1N matrice di iterazione
e(k)
e(k) 0e(0)P1N
7/23/2019 Metodi Analitici e Numerici Per l'Ingegneria Aerospaziale - Saleri
29/107
3. Risoluzione dei sistemi lineari
dove e il raggio spettrale
P(x(k+1) x(k)) = b+N x(k) P x(k)= b Ax(k)= r(k)
r(k) e il residuo, e una misura dellerrore che si sta commettendo
r(k) = 0 se x(k) =x
si puo pensare di introdurre un gdl, un acceleratore del metodo che potr aessere stazionario o dinamico k
P(x(k+1)
x(k)) = kr
(k)
es. se P = I
x(k+1) = x(k) +kr(k)
= x(k) +k(b Ax(k))= kb+ (I Ak)x(k)
P x = Nx+b
A = I NN = I Ax = (I
A)x+b
B = I kA (B)< 1si scelga k affinche (B) sia opportunamente piccolo
|(B)| = |1 k(A)|
Figura 3.2: Andamento del raggio spettrale.
0 < < 2
max
ott = 2
min+max
Appunti 27 F.Cadei BS
7/23/2019 Metodi Analitici e Numerici Per l'Ingegneria Aerospaziale - Saleri
30/107
3. Risoluzione dei sistemi lineari
3.2.1 Metodo del Gradiente
E equivalente a Richardson non stazionario, ma migliore. Lipotesi chesi avanza e quella di considerare solamente matriciA simmetriche definitepositive
Ax= b
equivale a trovare il minimo di un potenziale (x)
(x) = 1
2xTAx xTb
= 1
2[x1, . . . , xn]
a11 a1n...
. . . ...
an1 ann x1...
xn [x1, . . . , xn]b1...
bn =
1
2[x1, . . . , xn]
a11x1+ +a1nxn...an1x1+ +annxn
[x1b1+ +xnbn]=
1
2(a11x
21+a12x2x1+ +a1nxnx1+a21x1x2+ +annx2n)
(x1b1+ +xnbn)si sta cercando il minimo del funzionale, questa espressione va quindi derivata
rispetto ad xe posta uguale a zero
x1=
1
2(2a11x1+a12x2+ +a1nxn+a21x2+a31x3+ +an1xn) b1
ma A= AT aij =aji=
1
2(2a11+ 2a12x2+ + 2an1xn) b1
= a11+a12x2+ +an1xn b1 = 0e esattamente la prima equazione del sistema che si deve risolvere = 0 se x e soluzione del sistema; trovare il minimo del potenziale ecome risolvere il sistema lineare
partendo da un x(0) arbitrariox(k+1) =x(k) +d(k) k
uninterpretazione delle nuova posizione puo essere intesa come: partire da unpuntox(k), avanzare secondo una certa direzioned(k) per una certa quantitak; nel metodo di Richardson la direzione di spostamento era il residuo
(z) = Az bin x(k)
(x(k)) = Ax(k) b= r(k)
Appunti 28 F.Cadei BS
7/23/2019 Metodi Analitici e Numerici Per l'Ingegneria Aerospaziale - Saleri
31/107
3. Risoluzione dei sistemi lineari
ne risulta che il residuo e la direzione di massima discesa verso il minimo.Sulla direzione
x(k+1) = x(k) +kr(k)
x(k+1)
= 1
2x(k+1) TAx(k+1) x(k+1) Tb
= 1
2
x(k) +kr
(k)T
A
x(k) +kr(k) x(k) +kr(k)Tb
= 1
2x(k) TAx(k) +
1
2kr
(k) TAx(k) +1
2x(k) TAkr
(k) +
+1
2
kr(k) TAr(k)
x(k) Tb
kr
(k) Tb
A e simmetrica r(k) TAx(k) =x(k) TAr(k)=
1
2x(k) TAx(k) +kr
(k) TAx(k) +1
22kr
(k) TAr(k)
x(k) Tb kr(k) Tb= +kr(k) T(Ax(k) b
r(k)
) +
= 1
2x(k) TAx(k) kr(k) Tr(k) +1
22kr
(k) TAr(k) x(k) Tbsi deriva e azzera per ottenere lk che minimizza
ddk
= r(k) Tr(k) +kr(k) TAr(k) = 0
kott = r(k) Tr(k)
r(k) TAr(k)
si e ottenuta lespressione dinamica dik.Riassumendo: dato x(0) arbitrario, k 0 fino a convergenza si dovracalcolare
k=
r(k) Tr(k)
r(k) TAr(k)
x(k+1) =x(k) +kr(k)
r(k+1) =b Ax(k+1)
si calcolaksi aggiorna la soluzione
si aggiorna il residuo
Si puo pensare di ottimizzare il metodo, nel senso di ridurre il costo com-putazionale, gia da subito evitando di calcolare Ax(k+1) e riconducendo questocalcolo a conti gia eseguiti
r(k+1) = b A(x(k) kr(k))= b Ax(k) kAr(k)
r(k+1) = r(k) kAr(k)
Appunti 29 F.Cadei BS
7/23/2019 Metodi Analitici e Numerici Per l'Ingegneria Aerospaziale - Saleri
32/107
3. Risoluzione dei sistemi lineari
in questo modo invece di calcolare Ax(k+1) e Ar(k) si calcola solo la secondae la si utilizza piu volte.
Si tratti della convergenza: con la soluzione esatta il metodo si ferma
k = 0 r(k) = 0 x(k) =x
Th. di convergenzaIl metodo del Gradiente convergex(0) se A e simmetrica definita positiva einoltre
e(k+1)
2
K(A) 1K(A) + 1
k+1
e(0)
2
k
0
con: K(A) n di condizionamento della matrice, per A s.d.p. K(A) = maxmin
.Questo risulta essere un metodo che converge lentamente per una matricemal condizionata, ad esempio per K(A) = 100 labbattimento dellerrore eK(A)1K(A)+1
= 99101
, e vicino a 1, se fosse 1 si inchioda e non converge. Piu siimpiega a far convergere il metodo e piu gli errori hanno modo di propagarsi.
3.2.2 Metodo del Gradiente Coniugato
Il potenziale si puo pensare come una vallata di sezione trasversale piu
o meno ellittica e parabolica nella sezione longitudinale con al vertice lasoluzione esatta del problema di interesse. Per vallate lunghe e strettemax min
K(A) grande, paraboloide ellittico
K(A) piccolo, paraboloide tondo
per velocizzare il raggiungimento del minimo posso pensare di non ripeterele direzioni gia percorse durante la discesadirezioni ortogonali
r(k) Tr(k+1) = 0
r(j) Tr(k+1) ?= 0
si impone lortogonalita del residuo tra il passo k ed il passo k +1, quello checi si chiede e: tutte le direzioni risulteranno ortogonali? No.
r(k) Tr(k+1) =r(k) Tr(k) kr(k) TAr(k) = 0
localmente sono ottimale.
Appunti 30 F.Cadei BS
7/23/2019 Metodi Analitici e Numerici Per l'Ingegneria Aerospaziale - Saleri
33/107
3. Risoluzione dei sistemi lineari
Figura 3.3: Direzioni ortogonali.
r(j) Tr(k+1) = 0k non e tale da azzerare tutte le direzioni, risulteranno quindi ortogonali a2 a 2.
Figura 3.4: Direzioni ortogonali a coppie.
Si ricercano direzioni tra loro tutte ortogonali in quanto muoversi secondoqueste direzioni in uno spazio di dimensione n porta il metodo a convergenzain n passi.dato x(0) arbitrario, k
0 k=?
x(k+1) = x(k) +kd(k)
d(k) = ? dovranno essere, linearmente indip.Si cerchi ancora il minimo del potenziale
x(k+1)
= 1
2
x(k) +kd
(k)T
A
x(k) +kd(k) x(k) +kd(k)T b
= 1
2x(k) TAx(k) +kd
(k) T
Ax(k) b x(k) Tb+12
2kd(k) TAd(k)
d
dk x(k+1)= 0 k = d(k) Tr(k)
d(k) TAd(k)
d(k) si intenda usare, ricavato k resta da trovare d(k), si utilizza il residuoper trovare d in modo che siano ortogonali al residuo, il motivo e cher(k) euna combinazione lineare di d
r(k+1) = b Ax(k+1)= b Ax(k) +kd(k)= r(k) kAd(k)
Appunti 31 F.Cadei BS
7/23/2019 Metodi Analitici e Numerici Per l'Ingegneria Aerospaziale - Saleri
34/107
3. Risoluzione dei sistemi lineari
sed(k) Tr(k+1) =d(k) Tr(k) kd(k) TAd(k) = 0
la direzione d(k) risultaa r(k+1), se si riuscisse a trovare le nuove direzioniortogonali a tutte le precedenti si avrebbe la convergenza in n passi esatti.?d(j) Tr(k+1) = 0?, ci si chiede come devono essere fatte le d(j)
0 =d(j) Tr(k) kd(j) TAd(k)
si supponga di aver generatodral passo prima, allorad(j) Tr(k) = 0, affinchetutte le d siano ortogonali al passo seguente
d
(j) T
Ad
(k)
= 0
le direzioni devono essere ortogonali rispetto ad A, A-ortogonali oppure A-coniugate.Avviando la procedura: parto da una direzione ortogonale al primo residuo
d(0) = r(0) inf atti r(0) Tr(1) = 0 r(0)r(1)ora si imponga ched(1) sia A-ortogonale ad(0)
0 = d(1) TAd(0)
successivamente
d(2) t.c. d(2) TAd(0) = 0
d(2) TAd(1) = 0
d(k) ?
allaumentare dellindice continuano ad aumentare le condizioni da soddis-fare.SeA e simmetrica e definita positiva esiste un teorema che dimostra che percostruire direzioniA-ortogonali basta garantire le A-ortogonalita sullultimacoppia di direzioni perche si parte da direzioni furber(k+1), quello che si trovae A-ortogonale a tutto il resto
d(k+1) = r(k+1) kd(k)cerco k := d
(k) TAd(k+1) = 0
d(k) TAr(k+1) =kd(k) TAd(k)
k = d(k) TAr(k+1)
d(k) TAd(k)
k e un parametro che A-ortogonalizza rispetto allultima direzione, e unaprocedura che contiene gli errori di arrotondamento.
Appunti 32 F.Cadei BS
7/23/2019 Metodi Analitici e Numerici Per l'Ingegneria Aerospaziale - Saleri
35/107
3. Risoluzione dei sistemi lineari
Th. di convergenzaSe A e s.d.p. il metodo del Gradiente Coniugato converge alla soluzioneesatta in al piu n iterazioni3, con n =dim(A) e A Rnn.Il difetto e che questo metodo converge alla soluzione esatta in n passi sola-mente in aritmetica esatta; se si operano degli arrotondamenti literatannonha ancora spaziato tutto Rn; se si continua ad iterare con k > n si cominciaa fare danni. Il metodo si puo fermare impostando il numero massimo diiterate o controllando lerrore
e(k+1)A
K(A) 1
K(A) + 1
k+1e(0)A
e(k+1)A = e(k+1) TAe(k+1)si risente meno dellerrore grazie alla presenza della radice. Si puo pensare diridurre lerrore agendo su K(A) con tecniche di precondizionamento. Altricriteri darresto possono essere basati sul residuo, sulla prima iterata o sultermine noto
a. r(k) tolr(0)b. x(k+1) x(k)2 tole(0)
queste grandezze vanno valutate in un norma opportuna inoltre va consid-erato che rispetto agli altri metodi darresto la scalatura del residuo e piu
significativo perche nei termini di accelerazione k e k.
Figura 3.5: Andamento del numero di iterate in funzione di K(A).
3.3 Precondizionatori
I precondizionatori sono un modo per ridurre il numero di condizionamentodella matrice di partenza. Una rappresentazione geometrica del numero dicondizionamento e la seguente
3parte come metodo iterativo ed arriva come metodo diretto
Appunti 33 F.Cadei BS
7/23/2019 Metodi Analitici e Numerici Per l'Ingegneria Aerospaziale - Saleri
36/107
3. Risoluzione dei sistemi lineari
Figura 3.6: Rappresentazione geometrica di K(A).
considerando I ={insieme delle matrici singolari} si puo intendere il nu-mero di condizionamento come la distanza tra la matrice A e questo insieme;ora se consideriamo le approssimazioni introdotte dal sistema di calcolo (tron-camenti ecc.) la matrice A cambia avvicinandosi sempre piu allinsiemeI4.
Ax= b 1 sol
Ax= b sol0 sol
Si puo pensare di passare per via algebrica ad un problema equivalente (stessa
soluzione) ma meglio condizionato; partendo da
Ax= b
la strategia di pre-condizionamento trova P invertibile t.c.
K(P1A) K(A)
e risolve il sistemaP1Ax= P1b
TrovareP e un problema non lineare, la soluzione banale e data da P =A,
e inutile perche e vero che
K(P1A) = 1 K(A)
ma invertire P significa invertire Ache e il problema iniziale.
4qualora vi entrasse la A risulterebbe singolare, quindi non invertibile e non si puorisolvere il problema.
Appunti 34 F.Cadei BS
7/23/2019 Metodi Analitici e Numerici Per l'Ingegneria Aerospaziale - Saleri
37/107
3. Risoluzione dei sistemi lineari
Pnon deve essere solamente invertibile, deve anche essere facile da invertire
- Precondizionatore di Jacobi
P=PJ=diag(A)
se tutti gli ai= 0; e il piu semplice non funziona tranne quando lematrici sono a prevalenza diagonale stretta5
|aii| >n
j=1n=i
aij
- Precondizionatore in norma 2
P =diagAi2
con Ai si intende la riga i-esima di A; questa soluzione funziona parmatrici s.d.p. quando e forte la definita positivita
- Precondizionatore tri-diagonaleo di Gauss-Sider; funziona male ed inoltre si perde di simmetria
- Precondizionatore da fattorizzazione di Gauss incompleta
P A= LU rimpiazzoLU=Pcos invertendo Psi devono risolvere solo 2 sistemi triangolari6
- Precondizionatore polinomiale o di Neumann
A = P N=PI P1NA1 =
I P1N1P1
=
k=0 P1N
P1
troncando la serie si ottengono dei polinomi
= q
k=0
P1N
P1
metodo costoso converge per
P1N
< 1 ossia per q= 1, 2.
5per le differenze finite gia non va bene6si riduce loccupazione di memoria
Appunti 35 F.Cadei BS
7/23/2019 Metodi Analitici e Numerici Per l'Ingegneria Aerospaziale - Saleri
38/107
http:/
/aero.po
limi.tripod.c
om
Capitolo4Differenze Finite 2D e 3D
Il problema modello nel caso di problemi bidimensionali si presenta cometrovare u, R t.c.
u= 2ux2
2uy2
=f in
u= g su
deve essere rettangolo e f, g funzioni note1.Problema alle derivate parziali: si sostituisce con un insieme di punti, inodi di coordinate
(xi; yj )
xi = 0, 1, . . . , N
yj = 0, 1, . . . , M
Figura 4.1: Discretizzazione in nodi.
1problema di Laplace
36
7/23/2019 Metodi Analitici e Numerici Per l'Ingegneria Aerospaziale - Saleri
39/107
4. Differenze Finite 2D e 3D
u= f u(xi, yj) =f(xi, yj) nei pti interniu(xi, yj) =g(xi, yj) per (xi, yj)
lequazione funziona per i punti interni, la condizione al contorno per il bordo
2u
x2(xi, yj)
2u
y 2(xi, yj) =f(xi, yj)
separatamente si conoscono gia le approssimazioni
Figura 4.2: Numerazione a 2 indici.
2u
x2(xi, yj) = ui1,j 2ui,j+ ui+1,j
h2x
2u
y 2 (xi, yj) = ui,j1
2ui,j+ ui,j+1
h2y
ne segue
ui1,j 2ui,j+ui+1,jh2x
ui,j1 2ui,j+ui,j+1h2y
=f(xi, yj)
si cerca la sparsita della matrice 1
h2xui1,j 1h2y ui,j1+ 2ui,j
1h2x
+ 1h2y
1h2x
ui+1,j 1h2y ui,j+1=f(xi, yj)ui,j =g(xi, yj )
cambiando la numerazione si puo semplificare la struttura della matricecon questa tecnica si ottiene una matrice pentadiagonale del tiponel caso del problema completo a coefficienti costanti
u+x ux
+ yu
y+u = f
con u
x=
uI+1 uI12hx
con u
y =
uI+(N+1) uI(N+1)2hy
Appunti 37 F.Cadei BS
7/23/2019 Metodi Analitici e Numerici Per l'Ingegneria Aerospaziale - Saleri
40/107
4. Differenze Finite 2D e 3D
Figura 4.3: Numerazione a 1 indice e schema a croce.
Figura 4.4: Sparsita della matrice con lo schema a croce.
e se ci fossero derivate miste
2u
xy
(xI)
x
uy
(xI)
uy
(xI+1) uy (xI1)2hx
con u
y(xI+1)
uI+1+(N+1) uI+1(N+1)2hy
con u
y(xI1)
uI1+(N+1) uI1(N+1)2hy
si puo notare che le derivate miste allargano il supporto
Figura 4.5: Supporto per le derivate miste.
Appunti 38 F.Cadei BS
7/23/2019 Metodi Analitici e Numerici Per l'Ingegneria Aerospaziale - Saleri
41/107
4. Differenze Finite 2D e 3D
laccuratezza e 2 rispetto ad h con h = max(hx, hy); il passo variabile eimpiegato solitamente a blocchi e nei blocchih rimane costante.Per il problema
u = fw= uw= f
w1x
w2y
= f
w1 =u
xw2 =
u
y
si complica tutto.
Se e arbitrario e non rettangolare
Figura 4.6: Dominio generico.
il primo modo per discretizzare il problema e approssimare il bordo al reticolodi nodi che lo contiene interamente utilizzando la condizione al contorno sulnuovo bordo (che risultera squadrettato) e lequazione del problema sui nodiinterni. Nonostante lapprossimazione del bordo migliora affinando il passorisultano non definiti g ,,f (spesso valori sperimentali) fuori dal dominiofisico2.Si puo pensare di utilizzare un metodo di calcolo combinato che utilizzi dellegriglie non equispaziate per descrivere il bordo e griglie equispaziate per inodi interni.
2le condizioni al contorno sono valori sul bordo, il bordo squadrettato non e il bordodel problema
Appunti 39 F.Cadei BS
7/23/2019 Metodi Analitici e Numerici Per l'Ingegneria Aerospaziale - Saleri
42/107
4. Differenze Finite 2D e 3D
Figura 4.7: Soluzione alternativa.
Per questo tipo di approccio si pagano due codici: strutturato per i nodiinterni e non strutturato sul bordo, in compenso il metodo funziona.Unulteriore soluzione e la mappatura di
Figura 4.8: Mappatura del dominio.
si cerca una mappa che faccia lavorare sul rettangolo, perche sul dominiorettangolare le differenze finite sono ottimali. I problemi che si incontranosono legati allarrotolamento della discretizzazione, inoltre mappando si cam-biano i sistemi di riferimento e le derivate includono derivate incrociate e condirezione strane du
dx= dx
dxdudx +
dydy
dudy .
Sul vero dominio la griglia si deforma e puo creare accumulazioni inutili eaddirittura sforare il dominio stesso.
Figura 4.9: Deformazione del reticolo.
Appunti 40 F.Cadei BS
7/23/2019 Metodi Analitici e Numerici Per l'Ingegneria Aerospaziale - Saleri
43/107
4. Differenze Finite 2D e 3D
Nel caso di domini 3D tutto funziona come per i domini 2D, ma con unadirezione in piu
Figura 4.10: Discretizzazione del dominio 3D.
si suddivide in cubetti e la soluzione e solo sui vertici, anche in questi casise compaiono derivate miste il supporto dallo schema a diamante (7 nodi)diventa cubico.
Figura 4.11: Schema a diamante e diamante con derivate miste.
Anche per il 3D si puo pensare ad una tecnica di mappatura, uno degliimpieghi e lidrodinamica marina; si utilizzano le differenze finite mappateperche quando si ritorna al dominio fisico la mesh segue il fondo (Mappa ).
Figura 4.12: Mesh nel dominio fisico.
Appunti 41 F.Cadei BS
7/23/2019 Metodi Analitici e Numerici Per l'Ingegneria Aerospaziale - Saleri
44/107
http:/
/aero.po
limi.tripod.c
om
Capitolo5Spazi
5.1 Spazio L2
Definizione dello spazio L2 (e uno spazio di Lebesgue)
L2(0,1) =
v:
10
v2 dx <
Prodotto scalare in L2
(v, w)L2(0,1)
=
10
vw dx
Norma in L2
vL2(0,1)
=
(v, v)L2(0,1)
=
10
v2 dx
5.2 Spazi H1
Definizione dello spazio H1 (e uno spazio di Sobolev)
H1(0,1)=
v L2(0,1), v L2(0,1)
Definizione dello spazio H10
H10(0,1) =
v H1(0,1), v(0) =v(1) = 0
42
7/23/2019 Metodi Analitici e Numerici Per l'Ingegneria Aerospaziale - Saleri
45/107
5. Spazi
Prodotto scalare in H1 e H10
(v, w)H1(0,1) =
10
vw dx+
10
vw dx
Norma in H1 e H10
vH1(0,1)
=
(v, v)H1(0,1)
=
10
v2 dx+
10
v2 dx
u2H1 = u2L2+ u2L2
5.3 Spazi X1hDefinizione dello spazio X1h
X1h =
vh : vh|Ii P1, vh C0(0, L)
e un stottospazio di H1
Definizione dello spazio
X1h
X1h =
vh X1h, vh(0) =vh(L) = 0
5.4 Proprieta
Disuguaglianza di Cauchy-Schwartz
(v, w)L2(0,1) vL2(0,1) wL2(0,1)Disuguaglianza di Poincare
c >0 (funzione solo nellintervallo) t.c.u H10(0,L)uL2
(0,L) cuL2
(0,L)
Disuguaglianza tra spazi
uL2vL2 uVvVperche la norma in V (H1 o H10 ) contiene la stessa parte di quella in L
2 piuun altro pezzo positivo
Appunti 43 F.Cadei BS
7/23/2019 Metodi Analitici e Numerici Per l'Ingegneria Aerospaziale - Saleri
46/107
5. Spazi
5.5 Spazio di Hilbert
Caratteristiche dello spazio di Hilbert:
o Spazio vettoriale chiuso rispetto alla operazione somma, contiene lele-mento nullo, unita ecc.
o Spazio normato V R t.c.1.vV 0 evV = 0 v= 02. vV = || vV R,v V
3.v+wV vV + wV v, w Vo Il funzionale F : V R e lineare e continuoc > 0 t.c. |F(v)|
cvV v Vo Spazio in cui le successioni di Cauchy di elementi dello spazio conver-
gono a elementi dello spazio (spazi di Banach)
o Spazio chiuso rispetto al prodotto scalare (misuro gli angoli)
Appunti 44 F.Cadei BS
7/23/2019 Metodi Analitici e Numerici Per l'Ingegneria Aerospaziale - Saleri
47/107
http:/
/aero.po
limi.tripod.c
om
Capitolo6Formulazioni deboli
u =fmoltiplico per la funzione test e integro
1
0
uv dx=
10
fv dx
integro per parti
10
uv dx= 10
uv dx uv10
riscrivendo il problema 10
uv dx
uv1
0=
10
fv dx
Il problema debole risulta equivalente al problema variazionale:
cercare u
V :
J(u) =minvV
J(v) con
J(v) = 12 10(v)2 dx 10 fv dx
Si supponga che u sia soluzione del problema variazionale, allora ponendov= u+w si ha che
J(u) J(u+w) w Vsi ricerca il minimo, quindi mediante rapporto incrementale, si deriva emanda a zero
lim0
J(u+w) J(u)
= 0
45
7/23/2019 Metodi Analitici e Numerici Per l'Ingegneria Aerospaziale - Saleri
48/107
6. Formulazioni deboli
da cui
J(u+w) = 1
2
10
[(u+w)]2 dx 1
0
f(u+w) dx
= 1
2
10
[u2 +2w2 + 2uw]dx 1
0
fu dx 1
0
fw dx
= J(u) +1
2
10
[2w2 + 2uw]dx 1
0
fw dx
di conseguenza
J(u+w) J(u) =12 10 [w 2 + 2uw]dx 1
0fw dx
passando al limite e imponendo che si annulli si ottiene 10
uw dx 1
0
fw dx= 0
ovvero usoddisfa il problema debole.
6.1 Dirichlet Omogeneo
u(0) = u(1) = 0 restringo la ricerca alle v : v(0) = v(1) = 0 semplificandocon le condizioni al contorno il problema diventa
trovare u H10 t.c.v H10 10
uv dx=
10
fv dx
6.2 Dirichlet Non Omogeneo
u(0) =g0eu(1) =g1 sono le condizioni date al contorno riscrivo la soluzionenella forma
u(x) = u(x) +Rg
Riscrivendo il problema come soluzione del problema omogeneo
u =fu(0) = 0u(1) = 0
Appunti 46 F.Cadei BS
7/23/2019 Metodi Analitici e Numerici Per l'Ingegneria Aerospaziale - Saleri
49/107
6. Formulazioni deboli
e rilevamento del dato al contorno (definito perche e una retta)Rg = 0Rg(0) =g0
Rg(1) =g1
trovare u H10 t.c.v H10 1
0
u v dx=
10
fv dx 1
0
Rgv dx
6.3 Neumann Omogeneo
u+u= fformulazione del problema necessaria per avere una sola soluzioneu(0) = u(1) = 0 sono le condizioni date al contorno, moltiplicando per lafunzone test v e integrando per parti si ottiene 1
0
uv dx
uv1
0+
10
uv dx=
10
fv dx
semplificando con le condizioni al contorno il problema diventa
trovare u0 H1 t.c.v H1 10
uv dx+
10
uv dx=
10
fv dx
6.4 Neumann Non Omogeneo
u(0) = 0 e u(1) = 1 sono le condizioni date al contorno il probelma
assume la forma
trovare u H1 t.c.v H1 10
uv dx+
10
uv dx=
10
fv dx+1v(1) 0v(0)
Appunti 47 F.Cadei BS
7/23/2019 Metodi Analitici e Numerici Per l'Ingegneria Aerospaziale - Saleri
50/107
6. Formulazioni deboli
6.5 Funzionali e Forme
6.5.1 Funzionale
Il funzionale e una generalizzazione della funzione, e unapplicazione che vada uno spazio di funzioni a un numero reale
F :V Res.
F(v) = 1
2 1
0 |v
|2 dx
v
V
F(v) =
10
v dx
il primo esempio e un funzionale non lineare il secondo lineare, nel senso che
F(v+w) = F(v) +F(w)
, Rv, w V
Lintegrale ha ragion dessere in quanto distrugge la dipendenza da x restituen-do un numero.Riprendendo le formulazioni deboli si puo mostrare:
Dirichlet Omogeneo
F(v) =
10
fv dx
f L2(0,1) assegnata
v
H10 =V lineare
Dirichlet Non Omogeneo
F(v) =
10
fv dx 1
0
Rgv dx
f L2(0,1) assegnataRg H1(0,1)v H10 =V lineare
Appunti 48 F.Cadei BS
7/23/2019 Metodi Analitici e Numerici Per l'Ingegneria Aerospaziale - Saleri
51/107
6. Formulazioni deboli
Neumann Omogeneo
F(v) =
10
fv dx
f L2(0,1) assegnatav H1 =V lineare
Misura del Funzionale
Per valutare un funzionale gli si fa processare una funzione e se ne prende ilmassimo, un po come si e fatto per le matrici
|F| =sup |F(v)|||v||Vsi normalizza per svincolarsi dal caso particolare.
6.5.2 Forma
E unapplicazione binaria1 che va da una coppia di spazi ai numeri reali
a(, ) =V W R
ne e un esempio il prodotto scalare x
T
y R, si parte da due vettori e siarriva ad un numero.La differenza sostanziale tra Forme e Funzionali sta nel fatto che le formedanno lo stesso peso agli argomenti che le processano mentre nei funzionalila F e fissa e si da peso alla v 1
0
fv dx = a(f, v) f, v V 10
fv dx = F(v) v V
Riprendendo le formulazioni deboli
Dirichlet Omogeneo e Non Omogeneo
a(u, v) =
10
uv dx
u, v H10stessa forma ma cambia il funzionale
1tratta 2 argomenti
Appunti 49 F.Cadei BS
7/23/2019 Metodi Analitici e Numerici Per l'Ingegneria Aerospaziale - Saleri
52/107
6. Formulazioni deboli
Neumann Omogeneo e Non Omogeneo
a(u, v) =
10
uv dx+
10
uv dx
u, v H1
Forme bilineari
Per forme bilineare si intendono lineari rispetto entrambi gli argomenti
a(u+w, v) = a(u, v) +a(w, v)
, Ru,v,w Vlinearita I argomento
a(u,v+w) = a(u, v) +a(u, w)
, Ru,v,w Vlinearita II argomento
Forme simmetriche
a(u, v) =a(v, u) u, v V
6.6 Lemma di Lax-Milgram
Dato il problema
trovare u V t.c. a(u, v) =F(v)v Vla soluzione! se
1. V e uno spazio di Hilbert
2. la forma a: V V R e:o bilineare: a(u+w, v) =a(u, v) +a(w, v) e
a(u,v+w) =a(u, v) +a(u, w)
, R eu,v,w Vo continua:M >0 t.c.|a(u, v)| MuVvVv Vo coerciva: >0 t.c.|a(u, u)| u2Vu V
Appunti 50 F.Cadei BS
7/23/2019 Metodi Analitici e Numerici Per l'Ingegneria Aerospaziale - Saleri
53/107
6. Formulazioni deboli
3. il funzionale F :V
R e lineare e continuo:
c >0 t.c.
|F(v)| cvV v V
6.6.1 Verifica delle ipotesi
Continuita della forma a(u, v)
M >0 t.c. a(u, v) MuVvVa(u, v) =
L0
uv dx = (u, v)L2(0,L)
MuVvV
u
L2(0,L)v
L2(0,L) MuVvV
a(u, v) uVvV M= 1
Coercivita della forma a(u, v)
>0 t.c. a(u, u) u2Va(u, u) =
L0
(u)2 dx = u2L2(0,L)
conu2H1(0,L) = u2L2(0,L)
+ u2L2(0,L)
u
2
V C2
Pu
2
L2
(0,L)
+
u
2
L2
(0,L)
u2V (1 +C2P)u2L2(0,L) u2L2(0,L)
1
1 +C2Pu2V
a(u, u) 11 +C2P
u2V
Continuita del funzionale F(v)
c >0 t.c. F (v) cvV v VF(v) =
L0
fv dx = (f, v)L2(0,L) con f L2(0,L)
F(v) fL2(0,L) c
vL2(0,L)
u2V a(u, v) =F(v) cuV
Appunti 51 F.Cadei BS
7/23/2019 Metodi Analitici e Numerici Per l'Ingegneria Aerospaziale - Saleri
54/107
http:/
/aero.po
limi.tripod.c
om
Capitolo7Metodo degli Elementi Finiti
(Galerkin)
Trovareu V t.c. a(u, v) =F(v) v VIl problema viene risolto limitando lo spazio V delle funzioni test da spazioin infinite dimensioni a spazio finito di dimensione N(h) = dimVh < +,lapprossimazione e interna Vh
h0 V.Il problema che risolve Galerkin e
trovare uh Vh t.c. a(uh, vh) =F(vh) vh Vhapplicando il Lemma di Lax-Milgram ricavo la stabilita del sistema e la suacontrollabilita dai dati c e
! la soluzione uh e inoltreuh < CUtilizzo le basi dello spazio Vh per costruire u sugli elementi
uh(x) =
N(h)
i=1i(x)ui
a N(h)
i=1
uii(x), vh
= F(vh) vh Vh
N(h)i=1
ui a(i(x), vh) =F(vh) vh Vh
se vale per tutte le funzioni di base lho controllata per tutto lo spazio
vh V j B(Vh)
52
7/23/2019 Metodi Analitici e Numerici Per l'Ingegneria Aerospaziale - Saleri
55/107
7. Metodo degli Elementi Finiti (Galerkin)
N(h)i=1
ui a(i(x), j ) =F(j) j = 1, 2,...,N(h)
per qualsiasi Vh Galerkin mi rimanda ad un sistema lineare, basta che la
forma sia bilineare: Au = f conu =ui,
f =F(i), A= a(i, j)
Considero il problema modellou =f in (0, L); devo calcolare
a(i, j) =
L0
ij dx e F(j) =
L0
f j dx
prima scelgo le funzioni di base da utilizzare, utilizzo funzioni in X1h, polinomidi primo grado sullelemento con valore unitario sul vertice i-esimo e zero
altrove.i di base per X
1h, sul vertice sono
i(xj ) =
1 i= j
0 i =j
le funzionii(x) sullelemento risultano essere
i(x) =
xxi1
xixi1 se x (xi1, xi)xi+1xxi+1xi
se x (xi, xi+1)0 altrimenti
calcolo gli integrali
a(i, j) =
L0
ij =
0 se j (i 1, i , i+ 1)= 0 altrimenti
in questo modo A e tridiagonale, sarebbe piena se avessi usato le basi di
Lagrange.La forzante va approssimata e il funzionale calcolato con formule di quadratu-ra numerica;
F(i) =
L0
f i dx=
xi+1xi1
f i dx=ch
f(x) =
i
f(xi)i(x)
j
f(x)
xi+1xi1
j i dx
Appunti 53 F.Cadei BS
7/23/2019 Metodi Analitici e Numerici Per l'Ingegneria Aerospaziale - Saleri
56/107
7. Metodo degli Elementi Finiti (Galerkin)
Considero il caso di elementi tutti della stessa dimensione h
ai,i1 = a(i1, i) = xi
xi1
i1i dx=
xixi1
1h1
h dx= 1
h
ai,i+1 = ...= 1h
ai,i = a(i, i) =
xi+1xi1
(i)2 dx=
2
h
il sistema lineare che si ottiene e
A= 1
h
h 0 01 2 1
...... 1 2 1 ...
... . . . . . . . . .
...... 1 2 10 0 h
il vettore delle incognite e del termine noto sono
u=
u0u1
...
...uN(h)
uN(h)+1
f=
g0
f 1
...
...f N(h)
g1
porre attenzione al calcolo del termine noto, se uso la formula di quadraturadel trapezio ottengo
f i hf(x)
Lerrore di consistenza risulta:
h = a(u, vh)
F(vh)
= a(u, vh) a(uh, vh)= a(u uh, 0)
h = 0
il metodo e fortemente consistente perche non e necessario h 0 per avereerrore di consistenza nullo.Considerando la forma a(u uh, 0) = 0 come prodotto scalare, il lemma diCea rivela che u uh eVh sono ortogonali, il metodo di Galerkin e la migliorapprossimazione diu, u uh Vh.
Appunti 54 F.Cadei BS
7/23/2019 Metodi Analitici e Numerici Per l'Ingegneria Aerospaziale - Saleri
57/107
7. Metodo degli Elementi Finiti (Galerkin)
Figura 7.1: Significato geometrico del Lemma di Cea.
7.1 Convergenza
? u uh2V h0 0 ?
u uh2V a(u uh, u uh) coercivita della forma bilinearea(u uh, u vh+vh uh)a(u uh, u vh) +a(u uh, vh uh)a(u uh, u vh) + a(u uh, wh)
0 per consistenza forte
a(u uh, u vh)
u uh2
V Mu uhVu vhV continuita della forma bilineareu uhV M
u vhV vh Vh
se Vhh0 V allorau uh 0 per h 0
7.1.1 Ordine di convergenza
Le relazioni che reggono lordine di convergenza sono
u
r
h
u
L2
chl+1
u
Hs+1
u rhuH1 chluHs+1
con l = min(r, s), infatti utilizzando polinomi interpolanti di grado uno ecalcolando la norma in uno spazio di regolarita superiore a uno, non ci si pu oaspettare H1 di ordine maggiore al primo. Si nota inoltre che lerrore innorma L2 converge con un ordine in piu, cio e dovuto al fatto che H1 chiedeuna regolarita maggiore, per lappunto sulla derivata prima e derivando sidivide per h.
Appunti 55 F.Cadei BS
7/23/2019 Metodi Analitici e Numerici Per l'Ingegneria Aerospaziale - Saleri
58/107
7. Metodo degli Elementi Finiti (Galerkin)
7.2 Problemi di Diffusione-Trasporto(-Reazione)
u +u +u = f in (0, L)con coefficiente di diffusione o viscosita, coefficiente di trasporto e coefficiente di reazione; considero >0 e , qualsiasi.Considero il caso Dirichlet omogeneo e passo alla formulazione debole
L
0
uv dx+
L0
uv dx+
L0
uv dx =
L0
fv dx v H10 (0, L)
L
0
uv dx
uv
L
0 0 vH10+
L
0
uv dx+
L
0
uv dx =
L
0
fv dx
scelto lo spazio Vdevo verificare che tutti gli integrali esistano,uso Lemma di Lax-MilgramH10 e uno spazio di Hilbert
se f L2 (0, L) L0
fv dxcontinuita dia(u, v)
|a(u, v)|
L0
uv dx
+ ||
L0
uv dx
+
L0
uv dx
applico la disuguaglianza di Cauchy SchwartzuL2vL2+ ||uL2vL2+ uL2vL2
maggioro con norma in V e raccolgo
|a(u, v)| (+ || +) M
uVvV
coercivita di a(u, v)
|a(u, v)| u2V u V
L0
(u)2 dx+ || L
0
uu dx+
L0
u2 dx
u2L2 >0
+ 2 L
0
(u2) dx =0 Dirichlet omog.
+ u2L2 >0 trascuro
|a(u, v)| 1 +CP
u2V
per Lax-Milgram ! soluzione, so anche che
uVfL2
Appunti 56 F.Cadei BS
7/23/2019 Metodi Analitici e Numerici Per l'Ingegneria Aerospaziale - Saleri
59/107
7. Metodo degli Elementi Finiti (Galerkin)
dipende unicamente da , se
0 non ho piu il controllo sulla normadella soluzione, Galerkin funziona per diffusione dominante.Applico Galerkin
u uhV M(, , )()
u vh vh Vh
per M
grande non controllo lerrore, questo avviene per problemi a trasportoo reazione dominante:|| o .
7.3 Problema modello a trasporto dominanteu +u = 0u(0) = 1 u(1) = 0
cerco la soluzione fisica per = 1 e 0:se = 0 il problema diventa u = 0, u e una retta orizzontale, per > 0come in questo caso viene trascinata la soluzione al contorno sinistro per poiraccordarsi con la soluzione al contorno destro;se = riscrivendo il problema u + 1
u = 0 ottengo u = 0,u e ancora
una retta e congiunge le soluzioni al contorno.
Figura 7.2: Soluzione fisica del problema.
Se faccio processare il problema al calcolatore ottengo delle oscillazioni dellasoluzione che forniscono valori non fisici; questo accade fino a quandoh ,fastidioso.Applico Galerkin al problema per trovare quando non ho oscillazioni spurie,
Appunti 57 F.Cadei BS
7/23/2019 Metodi Analitici e Numerici Per l'Ingegneria Aerospaziale - Saleri
60/107
7. Metodo degli Elementi Finiti (Galerkin)
calcolo gli integrali per il termine di trasporto
uhi
ui1
i1i= 2 ui1
ui
ii= 0
ui+1
i+1i=
2ui+1
lo schema risulta
h(ui+1 2ui+ui1) + 2 (ui+1 ui1) = 0
u(0) = 1 u(L) = 0
raccolgo per termini simili
ui+1( h
+
2) + 2
hui+ui1(
h
2) = 0
Definisco il numero di Peclet (simile a Reynolds)
P e=||
2L P eh =
||h2
doveP eh e il numero di Peclet locale (sullelemento).
Moltiplicando la riformulazione del problema per h
ottengo
ui+1(P eh 1) + 2ui ui1(1 +P eh) = 0 i= 1, 2, . . . , N (h)sfruttando lanalogia con le equazioni alle differenze ottengo
ui = c1ri0+c2r
i1
simile c1e1t+c2e2t
= ri+1k (P eh 1) + 2rik ri1k (1 +P eh)= semplifico perr i1k= r2k(P eh 1) + 2rk (1 +P eh)= r0,1=1 P eh
P eh 1ui = c1+c2(
1 +P eh1 P eh )
i
ricavo le costanti con le condizioni al contorno, la soluzione esplode quandoesplode il termine ( 1+P eh
1P eh)N(h)+1, non ho oscillazioni spurie per P eh
7/23/2019 Metodi Analitici e Numerici Per l'Ingegneria Aerospaziale - Saleri
61/107
7. Metodo degli Elementi Finiti (Galerkin)
7.3.1 Decentramento del metodo di Galerkin (UpWind)
Voglio decentrare lo schema di Galerkin perche nei problemi a trasportodominante ha piu importanza quello che succede prima di quello che succededopo.Partendo dallo schema delle differenze finite decentrate indietro cerco unlegame con lo schema degli elementi finiti per decentrare indietro il problema
ui ui1h
= ui
h ui1
2h ui1
2h
= ui
h ui1
2h +
ui+12h
ui+12h
ui12h
= 12
2ui ui1 ui+1h
+ ui+1 ui12h
ui ui1h
= h2
ui+1 2ui+ui1h2
+ui+1 ui1
2h
il nuovo problema alle differenze finite diventa
(+ h2
)ui+1 2ui+ui1
h2 +
ui+1 ui12h
= 0
numerica = +
h
2
= (1 + h
2)
num. = (1 +P eh)
non produce oscillazioni spurie per qualsiasi e perP e >0 perche il numerodi Peclet del nuovo problema e sempre minore dellunita
P e
nuovo problema=
h
2num.
= h2(1 +P eh)
P e = P eh1 +P eh
7/23/2019 Metodi Analitici e Numerici Per l'Ingegneria Aerospaziale - Saleri
62/107
7. Metodo degli Elementi Finiti (Galerkin)
questo non e un metodo di Galerkin perche modifico il problema, inoltresarebbe valido solo nel caso di, costanti eh uniforme. Il metodo da unac-curatezza al piu di ordine uno; in realta sto risolvendo un altro problema, ilmio schema e come se vedesse un problema centrato con altra viscosita, lesoluzioni sono vicine a quelle del mio problema quando hsi avvicina allordinedi grandezza della viscosita fisica; in piu dimensioni butto via tutto.
7.4 Problema modello a reazione dominante
u +u = 0u(0) = 1 u(1) = 0
cerco la soluzione fisica per = 1 e 0:
Figura 7.3: Soluzione fisica del problema.
analogamente al caso di trasporto dominante valuto la soluzione per = 0 e= con le differenze finite ottengo la soluzione analitica, con gli elementifiniti non succede, perche la matrice di massa non e diagonale
uhvh = L
0
uhN(h)+1j=0 ujj
i dx
= (ui1
L0
iii dx+ui
L0
2i dx+ui+1
L0
i+ii dx)
= (ui1
xixi1
iii dx+ui
xi+1xi1
2i dx+ui+1
xi+1xi
i+ii dx)
= utilizzo la quadratura di Simpson
uhvh = (ui1
h
6+ui
2
3h+ui+1
h
6)
Appunti 60 F.Cadei BS
7/23/2019 Metodi Analitici e Numerici Per l'Ingegneria Aerospaziale - Saleri
63/107
7. Metodo degli Elementi Finiti (Galerkin)
e stata impiegata la quadratura di Simpson1 perche il trapezio da integralenullo.Lo schema risulta
h
(ui+1 2ui+ui1) + (ui1 h6
+ui2
3h+ui+1
h
6)
M matrice di massa
= 0
u(0) = 1 u(L) = 0
per = 0 trovo
ui1+ui+ui+1 = 0
u(0) = 1 u(L) = 0ho la soluzione a dente di sega ed e dovuto al fatto che la matrice di massanon e diagonale.
7.4.1 Condensazione della matrice di massa (Mass-Lumping)
Approssimo M per quadratura numerica in modo da renderla diagonale, laprocesso con funzioni a capanna fi X1h ne segue
mij = L
0
ij dx= = 0 se j (i 1, i , i+ 1)
= 0 altrimenti
mij =
0 se j=i= sej=i
in questo modo diagonalizzo M, inoltre volgio cheM M = o(h) ossiache lerrore vada almeno come lerrore di discretizzazione, uso la quadraturadel trapezio perche usa direttamente i nodi di discretizzazione
mij = ji12
(xi xi1)[i1(xi1)i(xi1) +i1(xi)i(xi)] = 0 j=i
xi
xi12i +
xi+1xi
2i = 12
h
2
1 =h = mii
12
(xi+1 xi)[i+1(xi+1)i(xi+1) +i+1(xi)i(xi)] = 0 j=isi verifica facilmente che
mii=1
j=1
mi,i+j
lelemento diagonale e la somma degli elementi della riga della matrice dimassa di partenza.
1f=
k
h
6[f(xk1) + 4f(xk1
2
) +f(xk)]
Appunti 61 F.Cadei BS
7/23/2019 Metodi Analitici e Numerici Per l'Ingegneria Aerospaziale - Saleri
64/107
7. Metodo degli Elementi Finiti (Galerkin)
Il discorso funziona anche nel caso in cui si utilizzino f
X2h, in questocaso la formula di quadratura piu appropriata e quella di Simpson in quantoquadro su 3 punti.Anche in questo caso esco dal metodo di Galerkin perche non uso la formaed il funzionale esatti (motivo per cui Galerkin in teoria non e praticabile).
Appunti 62 F.Cadei BS
7/23/2019 Metodi Analitici e Numerici Per l'Ingegneria Aerospaziale - Saleri
65/107
http:/
/aero.po
limi.tripod.c
om
Capitolo8Discretizzazione di ordine 2 in spazio
dipendenti dal tempo
Problema parabolico
u
t u
2ux2
+ uux
+u = f u= u(x, t)
x (0, L) e t [0, T]mi occupo di coefficienti variabili al massimo dipendenti dal solo spazio:
=(x) =(x) = (x) f=f(x, t)
il dominio computazionale e Q = (0, L) [0, T]Ho inoltre bisogno di condizioni iniziali e finali oltre che condizioni al bordo;per quanto riguarda le condizioni finali si lascera che il fenomeno venga guida-to dalla forzante f (problema di controllo) non fissando condizioni finali cherenderebbero il problema mal posto; quelle iniziali sono del tipoper t= 0 u(x, 0) =u(x);
per quanto riguarda le condizioni al bordo posso utilizzare condizioni di tipo:Dirichletu(x, t) =g(x, t), Neumann u(x, t) =(x, t) o miste.
8.1 Formulazione deboleu
t +Lu= f
e la formulazione forte, per passare alla debole, moltiplico per la funzionetest ed integro sul dominio Q
63
7/23/2019 Metodi Analitici e Numerici Per l'Ingegneria Aerospaziale - Saleri
66/107
8. Discretizzazione di ordine 2 in spazio dipendenti dal tempo
Tipo 1 Q
(u
t +Lu)v dQ=
Q
fv dQ
significa che anchev = v(x, t), se faccio cos sembra chet non abbia verso, lesoluzioni dipendono dal anche dal futuro.Tipo 2fisso il tempo T
t [0, T]cerco u(x, t) V t.c.v= v(x) V
L0
u
t v dx+ L0 Luv dx= L
0fv dx
affetto Q per ogni istante temporale cos non integro su t.A tempo congelato conosco tutto L
0
u
tv dx+a(u, v) =F(v)
V? non dipende da t
V =
H1 (0, L)
H10 (0, L)
Ho soluzione quandot (0, T] vale lhp di Lax-Milgram ! della soluzionela coercivita diventa
ut
v u=v 1
2t
u2 = 1
2t
u2L2 u(t, 0) = u0(x) indebolisco la coercivita, hotralasciato t= 0 quindi fisso una condizione al contornox (0, L)coercivita debole:
a(u, u) +u2L2 u2V
8.2 Discretizzazione in spazio con ElementiFiniti
Galerkin elementi finitit [0, T], devo dare una dimensione finita a V
Vh t.c. dimVh = N(h)< + e Vh h0 V
ne segue che Vh e uno spazio X1h
X1h
Xrhmi gioco lo spostamento dei nodi nel tempo, ma io voglio V fisso, in questo
Appunti 64 F.Cadei BS
7/23/2019 Metodi Analitici e Numerici Per l'Ingegneria Aerospaziale - Saleri
67/107
8. Discretizzazione di ordine 2 in spazio dipendenti dal tempo
modo si ottengono matrici di massa e rigidezza indipendenti dal tempo e sipossono calcolare una volta sola.
t (0, T]cerco uh=uh(x, t) Vh t.c. L0
uht
vh dx+a(uh, vh) =F(v) vh Vhuh(x, 0) =u0(x) x (0, L)
in generale non posso scrivere u0(x) in quanto e una funzione continua, vaquindi discretizzata, meglio scrivere uh(x, 0) = hu0(x); e va come lerroredellapprossimazione del problema
uh(x, t) =
N(h)i=1
ui(t)i(x)
ui(t) perche devo sentire il tempo in qualche modo,i(x) percheVh=V(x, t)
vh(x) =i(x) t (0, T) trovare per i= 1,...,N(h) L0
t
N(h)j=1
uj(t)j(x)
i(x) dx+aN(h)
j=1
uj(t)j(x), i
=
L0
f(x, t)i(x) dx
N(h)j=1
ujt
L0
j (x)i(x) mij
dx+N(h)j=1
uj(t) a(j(x), i) kij
= L
0
f(x, t)i(x) dx
uj = uj(t) dipende solo dal tempo, riscrivendo u(t) = uj(t) il problemadiventa
Mdudt
+Ku= F(t) t (0, T]u(0) =u0,h
dove u0,h e linterpolato i-esimo u0,ni = hu0(ui), si ottiene un sistema di
equazioni differenziali ordinarie, con M e Kche non dipendono dal tempo,vengono calcolate una volta sola per aver congelato la griglia di calcolo.
Appunti 65 F.Cadei BS
7/23/2019 Metodi Analitici e Numerici Per l'Ingegneria Aerospaziale - Saleri
68/107
8. Discretizzazione di ordine 2 in spazio dipendenti dal tempo
8.3 Discretizzazione in tempo con DifferenzeFinite
8.3.1 Eulero Esplicito
esplicito perche dipende solo dallistante prima (derivo indietro), il problemadiventa
u(tn) un t.c.Mu
n1un
t +Kun = F(tn) n 0
u0 = u0,h
risolvendoMun+1 =M un tKun + t F(tn)
ad ogni passo devo risolvere un sistema lineare, M non dipende dal tempofattorizzo una volta per tutte e risolvo ad ogni passo sistemi triangolariM=RTR Cholesky perche M e sdpoppure sostituiscoM con M usando mass lumping, riduco la dimensionedellinversione allinversione di una matrice diagonale.
8.3.2 Eulero Implicitoimplicito perche dipende, oltre che dallistante precedente, anche da quelloattuale, il sistema e non lineare (derivo in avanti)
Mun+1un
t +Kun+1 = F(tn+1)
u0=u0,h
1t
M+K non e diagonale, ne diagonalizzabile, posso solo fattorizzare conCholesky o Gauss, quando le matrici sono grandi uso sistemi iterativi. E.Eed E.I sono accurati di ordine 1, lerrore totale del sistema spazio-tempo e
una combinazione lineare dei due errori: ertot=chr + ct la parte riguardanteil tempo e di ordine 1 mentre r puo aumentare,
8.3.3 Crank-Nicolson
cerco lordine 2, combino E.E ed E.I. e genero i metodi[0, 1]
Mun+1 un
t +Kun+1 + (1 )Kun = F(tn+1) + (1 ) F(tn)
Appunti 66 F.Cadei BS
7/23/2019 Metodi Analitici e Numerici Per l'Ingegneria Aerospaziale - Saleri
69/107
8. Discretizzazione di ordine 2 in spazio dipendenti dal tempo
= 1 Eulero Implicito
= 0 Eulero Esplicito
= 1/2 Crank Nicolsonoss.1 tutti i metodi della famiglia tranne per = 0 sono impliciti.oss.2 tutti convergenti di oridne 1 eccetto per = 1/2 che e di ordine 2.per Crank-Nicolson ho un difetto, ho delle oscillazione per dati iniziali dis-continui, per questo motivo non uso esattametne = 1/2 quindi avrei ordine1, ma ho una costante di convergenza alta.
Figura 8.1: Andamento dellerore rispetto alla discretizzazione.
8.3.4 Analisi di stabilita per il -metodo
Pb. di Cauchy
y(t) =f(t, y(t)) t >0
y(0) =y0
Per Zero stabilita si intende stabilita su intervalli limitati, ossia in unintervallo limitato la soluzione non deve andare allinfinito.
|y(t)
y(t)
| c
|y0
y0
|Con lAssoluta stabilita si pretende che
t0>0 : t 0
y(0) =y0
Appunti 67 F.Cadei BS
7/23/2019 Metodi Analitici e Numerici Per l'Ingegneria Aerospaziale - Saleri
70/107
8. Discretizzazione di ordine 2 in spazio dipendenti dal tempo
C quindi se R()< 0 y(t) =et e
limt
y(t) = 0
Assoluta Stabilita per Eulero Implicito
1
t(un+1h unh, vh) +a(un+1h , vh) = (fn+1, vh) vh Vh
unh =L
0(un+1h unh)vh
prendendo f 0 elimino la parte non linearevh=u
n+1
h
e la soluzione al tempo attuale
1
t(un+1h , u
n+1h ) +a(u
n+1h , u
n+1h ) =
1
t(unh, u
n+1h )
per disug. Cauchy Schwartz1
tun+1h 2L2(0,L)+a(un+1h , un+1h )
1
tun+1h L2(0,L)unhL2(0,L)
per disug. Y ang [2ab a2 +b2] 1
t
1
2
un+1h 2L2(0,L)+ unh2L2(0,L)per la coercivita
1t
un+1h 2L2(0,L)+un+1h 2V per def. di V
1t
un+1h 2L2(0,L)+un+1h 2L2(0,L) riformulando
1
2tun+1h 2L2(0,L)+un+1h 2L2(0,L)
1
2tunh2L2(0,L)
1 un+1h 2L2(0,L) unh2L2(0,L) come per la soluzione analitica per f = 0 lasoluzione e controllata dai dati iniziali
2 un+1h 2L2(0,L) 11 + 2t
7/23/2019 Metodi Analitici e Numerici Per l'Ingegneria Aerospaziale - Saleri
71/107
8. Discretizzazione di ordine 2 in spazio dipendenti dal tempo
se
, w e lautofunzione rispetto allautovalore cercoxt.c. yTAx= yTx y; hoPasso al discreto
a(wh, vh) K rigidezza
=h(wh, vh) vh Vh
ih ne ho in numero finito = alla dimensione dello spazio
K wh= hM wh
se inverto Mho la ricerca di autovalori classica
(M1K)wh = hwh
in questo modo ho tanti h quanto la dimensione dello spaziovoglio a( , ):simmetrica, per godere delle proprieta della matrice A per lAsintotica Sta-bilita;
a(v, u) =a(u, v) u, v Vlinearmente indipendenti{ wih}N(h)i , forniscono una base per Vh= Vh = xiche e complicata da calcolareortogonale,
(wjh, wih)L2(0,L)=ji =
1 se j =i
0 se j=i(wjh, w
ih) nel senso degli Elementi Finiti, M=I
oss.1
un+1h =
N(h)j=i
un+1j wjh(x)
1t
(
N(h)j=i
un+1j wjh(x), w
ih(x)) +
a(N(h)j=i
un+1j wjh(x), w
ih(x)) +
(1 ) a(N(h)j=i
unj wjh(x), w
ih(x)) =
1
t(
N(h)j=i
unj wjh(x), w
ih(x))
Appunti 69 F.Cadei BS
7/23/2019 Metodi Analitici e Numerici Per l'Ingegneria Aerospaziale - Saleri
72/107
8. Discretizzazione di ordine 2 in spazio dipendenti dal tempo
ho il-metodo con sostituite le funzioni a capanna con quelle nuove= 0 solo per i= j , semplificazioni dovute allortonormalita
a(wj, wi) =j(wjh, w
ih)
1
tun+1i (w
jh, w
ih)
1
+ un+1i i+ (1 ) unii= 1
tuni
riordinandoun+1h (1 + ti) =u
ni(1 (1 )ti)
un+1h =1 (1 )ti
1 + tiunh i= 1 . . . N (h)
oss.2 voglio un+1h 0 per n 1 (1 )ti1 + ti
>0
< 1|1 (1 )ti|
7/23/2019 Metodi Analitici e Numerici Per l'Ingegneria Aerospaziale - Saleri
73/107
8. Discretizzazione di ordine 2 in spazio dipendenti dal tempo
si possono aggiungere altre perturbazioni durante la risoluzione (arrotonda-menti ecc.).sono zero-stabile se riesco a controllare la soluzione a partire dai dati, se:
|0|, |1|, . . . , |n| |un zn| c
volere che si verifichi con tutti i t e esageratomi basta che si verifichi
t t0zero-stabilita: proprieta di essere controllato dai dati; es. il valore deglizeri delle parabole, mentre il numero di radici non e controllabile dai daticontinuamente, ad un certo punto scatto da 2 a 1 a 0.I metodi a un passo sono zero-stabili, se:
- f C0 rispetto ad entrambi gli argomenti- f
y C0
oppure
- fLipschitziana rispetto al 2 argomento sec >0|f(t, y(t)) f(t, z(t))| c|y(t) z(t)|due derivate finite, e Lipschitziana ma non derivabile.
sottraggo i primi elementi dei sistemi
|un+1 zn+1| = |un zn| + t[f(tn, un) f(tn, zn) c|unzn|
]
|un+1 zn+1| = (1 +ct c
)|un zn|
Lasintotica stabilita e la stabilita su intervalli illimitati
y
(t) =f(t, y(t)) y(t)
t >0
y(0) =y0
y(t) e la scelta del problema a cui applicare la soluzione con C R()< 0(non e un esempio, e il problema sul quale si studia lassoluta stabilita)
y(t) =y(t)
y(0) =y0y(t) =et
t 0
Appunti 71 F.Cadei BS
7/23/2019 Metodi Analitici e Numerici Per l'Ingegneria Aerospaziale - Saleri
74/107
8. Discretizzazione di ordine 2 in spazio dipendenti dal tempo
Figura 8.2: Scelta del problema per lo studio della asintotica stabilita.
metodo numerico e asintoticamente stabile se (riproduce allinfinito la soluzione)t t0 accade che lim
nun= 0
Eulero Esplicito e asintoticamente stabile
un+1 = un+hf(tn, un)
un+1 = un+hun
= un(1 +h)
= u0(1 + h
z)n+1
Figura 8.3: Zona di asintotica stabilita.
in questo caso h= t per
0< t < 2
||condizione di asintotica stabilita per Eulero Esplicito
Appunti 72 F.Cadei BS
7/23/2019 Metodi Analitici e Numerici Per l'Ingegneria Aerospaziale - Saleri
75/107
8. Discretizzazione di ordine 2 in spazio dipendenti dal tempo
Eulero Implicito
un+1 = un+ tun+1
= un1 t
= u0
(1 t)n+1sempre< 1, incondizionatamente asintoticamente stabile tperche R() 0 e si dica se il metodo del gradienteconiugato puo essere applicato per la risoluz
Top Related