Laurea Magistrale in Ingegneria per l’Innovazione del Prodotto

Corso di Calcolo Numerico

Il metodo degli Elementi Finiti

Indice

1 Spazi funzionali lineari 1

2 Metodi variazionali 4

2.1 Il metodo di Galerkin . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4

2.2 Formulazioni deboli . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 7

3 Il metodo degli elementi finiti 9

3.1 Elementi finiti mono-dimensionali . . . . . . . . . . . . . . . . . . . . . . . . . . . . 9

3.2 Elementi finiti lagrangiani . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 13

3.3 Elementi finiti triangolari . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 17

4 Soluzione agli elementi finiti dell’equazione di Poisson 19

4.1 Caso 1-D . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 19

4.2 Caso 2-D . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 23

Esercizi 26

I

1 Spazi funzionali lineari

Uno spazio lineare e un’entita matematica che descrive un insieme di oggetti dotati di particolari

proprieta. In particolare, detto X = x, y, z, . . . tale insieme, esso definisce uno spazio lineare

reale se:

1. per ogni x e y appartenenti a X , e possibile definire l’elemento somma x + y anch’esso

appartenente a X ;

2. per ogni scalare α ∈ R, e possibile definire il prodotto αx appartenente a X ;

3. le operazioni precedentemente definite godono della proprieta commutativa, associativa e

distributiva;

4. per ogni x e y appartenenti a X , e possibile definire un prodotto interno il cui risultato sia

uno scalare in R.

Uno spazio lineare reale e, ad esempio, lo spazio vettoriale Rn i cui elementi sono costituiti da

ennuple di numeri reali ed il prodotto interno e il prodotto scalare. Se gli elementi di uno spazio

lineare reale sono funzioni, allora si parla di spazio funzionale lineare. In questo senso, uno spazio

funzionale puo essere interpretato come una generalizzazione del concetto di spazio vettoriale con

analoghe proprieta.

Quando una combinazione lineare di n elementi di uno spazio lineare:

α1x1 + α2x2 + . . .+ αnxn (1)

da l’elemento nullo se e solo se α1 = α2 = . . . = αn = 0, allora gli elementi x1, x2, . . . , xn

si dicono linearmente indipendenti. Qualora qualsiasi combinazione lineare di n + 1 elementi

appartenenti a X sia dipendente, si dice che la dimensione dello spazio e n e che ogni insieme di

n elementi linearmente indipendenti ne costituisce una base. Viceversa, se per ogni n > 0 esiste

sempre un insieme di n elementi linearmente indipendenti, si dice che lo spazio X ha dimensione

infinita. Gli spazi vettoriali sono esempi di spazi a dimensione finita, mentre gli spazi funzionali

hanno dimensione infinita. Tale proprieta consente di generalizzare facilmente alcune operazioni

definite in uno spazio vettoriale interpretando una funzione come un vettore ad infinite componenti.

Cio implica, in sostanza, il passaggio da operazioni di sommatoria discreta ad integrazioni nel

continuo. Siano, ad esempio, u e v due vettori di Rn. Come noto, il loro prodotto scalare uTv si

definisce come somma dei prodotti delle componenti, cioe :

uTv =

n∑

i=1

uivi (2)

Due elementi u(x) e v(x) definiti in un dominio R ed appartenenti ad uno spazio funzionale lineare,

invece, sono dotati di infinite componenti e la definizione di prodotto scalare si puo facilmente

1

estendere dalla (2) come:

< u, v >=

∫

R

u(x)v(x) dx (3)

In modo analogo e possibile definire una nozione di norma intesa come misura della grandezza di

una funzione. Ad esempio, come in Rn la norma euclidea di un vettore u si definisce come:

‖u‖2 =(

uTu)1/2

(4)

cosı in uno spazio funzionale la norma euclidea di u(x) risulta:

‖u‖2 =

[∫

R

u2(x) dx

]1/2

(5)

Uno spazio funzionale in cui sia possibile definire una norma come in (5) e detto normato. Inoltre,

quando il prodotto scalare (3) e finito per ogni coppia di funzioni, lo spazio funzionale e detto

misurabile.

Dati due sottoinsiemi S1 e S2 di uno spazio lineare S, e possibile definire una trasformazione che

colleghi ogni elemento del primo insieme, detto dominio, ad un solo elemento del secondo insieme,

detto codominio. Chiamiamo operatore la rappresentazione simbolica di tale trasformazione:

Ax = y (6)

dove x ∈ S1, y ∈ S2 e A denota la regola che trasforma x in y. L’operatore A e detto lineare se

per ogni α e β reali vale che:

A (αx+ βy) = αAx+ βAy (7)

Ad esempio, nello spazio vettoriale Rn l’operatore A puo essere rappresentato da una matrice

quadrata di dimensione n con la trasformazione da esso attuata definita da un prodotto matrice-

vettore. In maniera analoga in uno spazio funzionale si possono definire degli operatori lineari

di tipo differenziale, contenenti cioe operazioni di derivazione parziale. Ad esempio, l’operatore

differenziale:

A =∂2

∂x2+∂2

∂y2(8)

e normalmente indicato come operatore laplaciano ed esegue una trasformazione lineare su una

funzione u(x, y) in f(x, y). L’equazione differenziale di Poisson:

∂2u

∂x2+∂2u

∂y2= f (9)

si puo pertanto interpretare come l’applicazione dell’operatore laplaciano ad u per ottenere f . Il

laplaciano in (8) e frequentemente indicato anche con il simbolo ∇2. L’equazione differenziale (9),

pertanto, si scrive anche come:

∇2u = f (10)

2

Data l’analogia con le matrici negli spazi vettoriali, un operatore differenziale lineare puo essere

simmetrico e definito positivo. Come una matrice A e simmetrica se per ogni u e v non nullo in

Rn vale:

uTAv = vTAu (11)

cosı, ricordando la definizione (3), un operatore differenziale A e simmetrico se per ogni funzione

u(x) e v(x) non nulla di un determinato spazio funzionale F vale:

< u,Av >=< v,Au > (12)

Analogamente, come una matrice A e definita positiva se per ogni u non nullo in Rn vale:

uTAu > 0 (13)

cosı un operatore differenziale A e definito positivo se per ogni funzione u(x) non nulla in F vale:

< u,Au >> 0 (14)

Si consideri ora una successione di funzioni ξ1, ξ2, . . . , ξn linearmente indipendenti appartenente

allo spazio funzionale misurabile e normato F . Poiche F ha dimensione infinita, nessun insieme

finito di funzioni indipendenti puo a rigore costituirne una base. E, tuttavia, desiderabile cercare di

estendere il concetto di base anche per gli spazi funzionali. Diciamo che la successione di funzioni

ξ1, ξ2, . . . , ξn e completa o chiusa se per ogni u(x) ∈ F e per ogni scalare ǫ > 0 esiste un numero

n e un insieme di reali α1, α2, . . . , αn tali che per la funzione:

un = α1ξ1 + α2ξ2 + . . .+ αnξn (15)

valga:

‖u− un‖ < ǫ (16)

dove il simbolo ‖·‖ indica la norma definita in F . In altri termini, possiamo dire che la successione

ξ1, ξ2, . . . , ξn definisce un sottospazio Fn di F di dimensione finita e che essa si dice completa o

chiusa se qualsiasi funzione di F puo essere approssimata sufficientemente bene in Fn. In tal caso

si dice che la successione ξ1, ξ2, . . . , ξn costituisce una base anche per F ed ogni ξi viene anche

chiamata funzione base o funzione forma. Si osservi che la (16) garantisce la convergenza in media

di un ad u, ma non la convergenza uniforme. Cio significa che un potrebbe essere significativamente

diversa da u in qualche punto isolato, pur rispettando la condizione (16).

A titolo di esempio, si consideri lo spazio funzionale definito da tutte le funzioni continue

e illimitatamente derivabili nell’intervallo [0, 1]. La successione di funzioni 1, x, x2, x3, . . . , xn

e completa in questo spazio e ne costituisce una base. Infatti, ogni funzione continua e illimitata-

mente derivabile in un intervallo chiuso e limitato puo essere sempre approssimata sufficientemente

3

bene da un polinomio. E intuitivo che la rappresentativita del sottospazio generato da queste fun-

zioni aumenta al crescere di n: tanto maggiore e n, tanto migliore sara il grado di approssimazione

di un rispetto ad u. Invece la successione 1, x2, x4, . . . , x2n non e completa perche non consente di

rappresentare le funzioni dispari nel dominio [0, 1].

2 Metodi variazionali

Si consideri il seguente problema differenziale lineare:

Au = f su R (17)

Gu = q su ∂R (18)

dove u ed f sono funzioni definite nel dominio R, q e definita lungo la frontiera ∂R, e A e G sono

degli operatori differenziali lineari. Il problema ai limiti costituito dalle (17) e (18) descrive nel

formalismo degli spazi funzionali la ricerca della soluzione u di un’equazione differenziale definita

in R con determinate condizioni al contorno. Risultera pertanto G = 1 lungo la porzione di

frontiera ∂R1 in cui si impongono condizioni principali (o di Dirichlet) e G = ∂/∂n lungo ∂R2

in cui si impongono condizioni naturali (o di Neumann). Se q = 0 la (18) definisce condizioni al

contorno omogenee.

Sia F lo spazio funzionale lineare che contiene tutte le funzioni definite in R che soddisfano

le condizioni al contorno (18). La soluzione u del problema (17)-(18) appartiene ad F . I metodi

variazionali cercano una approssimazione un di u in un sottospazio a dimensione finita Fn di F

come combinazione lineare della successione di funzioni inidpendenti ξ1, ξ2, . . . , ξn che ne costi-

tuiscono la base. Se tale successione e completa, la soluzione esatta u potra essere approssimata

sufficientemente bene da un a patto di scegliere un valore n abbastanza elevato. Poiche le funzioni

base ξ1, ξ2, . . . , ξn sono scelte dall’operatore, la soluzione approssimata del problema (17)-(18) si

riduce alla ricerca dei coefficienti α1, α2, . . . , αn della combinazione lineare (15) che minimizzano

in qualche modo la distanza di un da u.

2.1 Il metodo di Galerkin

Sia un definita dall’equazione (15) la soluzione approssimata al problema (17)-(18). In generale,

la funzione residuo rn:

rn = Aun − f (19)

risultera essere non nulla in R. Le n incognite α1, α2, . . . , αn devono essere determinate mediante

n condizioni tali da minimizzare in qualche modo rn in R.

Per introdurre la tecnica di minimizzazione adottata e utile tornare all’analogia con gli spazi

vettoriali. Si consideri per semplicita un vettore r definito in R2. Un modo per imporre che r

4

x

x

v

v

1

2

2

1

Figura 1: Esempio di due vettori indipendenti in R2.

abbia lunghezza minima, cioe sia nullo, e che, dati due vettori v1 e v2 indipendenti in R2 (Figura

1), valga:

rTv1 = 0

rTv2 = 0

(20)

vale a dire che r sia simultaneamente ortogonale a v1 e v2. E, infatti, immediato osservare che

l’unico vettore di R2 ortogonale sia a v1 che a v2 e il vettore nullo. Si noti che la condizione di

indipendenza di v1 e v2 e fondamentale. Se per assurdo v1 e v2 fossero dipendenti, cioe paralleli

in R2, esisterebbe un’infinita di vettori r non nulli contemporanemente ortogonali ad entrambi.

La condizione (20) evidenziata in R2 per via geometrica puo essere facilmente estesa a Rn.

In questo caso r e minimo, cioe nullo, se simultaneamente ortogonale a n vettori linearmente

indipendenti v1, v2, . . . , vn:

rTvi = 0 i = 1, . . . , n (21)

La minimizzazione di r in Rn mediante le equazioni (21) puo facilmente essere estesa al caso di

uno spazio funzionale lineare data l’analogia evidenziata nel precedente paragrafo. Se rn e una

funzione in Fn e v1, v2, . . . , vn una successione di funzioni linearmente indipendenti, la condizione

di minimo per rn viene ottenuta imponendo che:

< rn, vi >= 0 i = 1, . . . , n (22)

cioe, ricordando le definizioni (3) e (19):

∫

R

(Aun − f) vi dR = 0 i = 1, . . . , n (23)

Le n equazioni (23) costituiscono un sistema determinato nelle n incognite α1, α2, . . . , αn che, risol-

to, consente di calcolare la soluzione approssimata cercata. Le condizioni di ortogonalita imposte

5

mediante le (23) sono note come condizioni di Petrov-Galerkin. Poiche le funzioni v1, v2, . . . , vn co-

stituiscono una base per un nuovo sottospazio funzionale Ln, si dice anche che un viene determinata

in Fn ortogonalizzando il residuo a Ln.

Una delle scelte piu semplici per la successione di funzioni v1, v2, . . . , vn e rappresentata dal-

l’assumere:

vi = ξi i = 1, . . . , n (24)

dato che la successione ξ1, ξ2, . . . , ξn costituisce proprio un insieme di funzioni linearmente indi-

pendenti. In questo caso gli integrali (23) diventano:

∫

R

(Aun − f) ξi dR = 0 i = 1, . . . , n (25)

e prendono il nome di integrali variazionali di Galerkin. Il metodo variazionale che discende

dall’applicazione delle (25) e noto come metodo di Galerkin.

Nel caso in cui l’operatore differenziale A sia simmetrico e definito positivo secondo le definizioni

(12) e (14), la soluzione approssimata un soddisfa un’interessante proprieta ottimale. Definita la

funzione errore en come:

en = un − u (26)

la soluzione un e tale da rendere minima la quantita :

< en, Aen >=

∫

R

A (un − u) (un − u) dR (27)

nota come norma energia associata all’operatore A. Infatti, ricordando la (17) e la linearita di A,

l’integrale in (27) puo essere riscritto come:

Ω (un) =

∫

R

(Aun − f) (un − u) dR (28)

Poiche un, fissata la base ξ1, ξ2, . . . , ξn, dipende solamente dagli scalari α1, α2, . . . , αn utilizzati

nella combinazione lineare (15), per minimizzare l’integrale in (28) e sufficiente annullare tutte le

derivate di Ω fatte rispetto ad αj (j = 1, . . . , n):

∂Ω

∂αj=

∫

R

[Aξj (un − u) + (Aun − f) ξj ] dR = 0 j = 1, . . . , n (29)

Sfruttando ora la simmetria di A le equazioni (29) diventano:

∫

R

2 (Aun − f) ξj dR = 0 j = 1, . . . , n (30)

che coincidono con le (25). Pertanto se A e simmetrico e definito positivo la soluzione approssimata

un rappresenta anche il punto di minimo in Fn per l’energia associata ad A.

6

2.2 Formulazioni deboli

Gli integrali di Galerkin (25) possono essere riscritti come segue:

∫

R

(Aun) ξi dR =

∫

R

fξi dR i = 1, . . . , n (31)

dove il membro a destra del segno di uguaglianza, non essendo funzione di un, costituisce un termi-

ne noto. Evidenziamo la dipendenza dalle incognite nell’integrale a primo membro introducendo

in (31) la definizione (15) di un e sfruttando la linearita di A:

n∑

j=1

[∫

R

(Aξj) ξi dR

]

αj =

∫

R

fξi dR i = 1, . . . , n (32)

Le equazioni (32) mostrano che le funzioni base ξi devono essere derivabili un numero di volte

almeno pari all’ordine dell’operatore differenzialeA. Ad esempio, se A = ∇2 le funzioni ξi dovranno

essere derivabili almeno due volte rispetto ad x e y affinche l’integrale a primo membro nelle (32) sia

non nullo. Se, tuttavia, si applica opportunamente la formula di integrazione per parti e possibile

ridurre l’ordine massimo di derivabilita richiesto per le funzioni base.

Mostriamo questo risultato con un esempio in un problema differenziale mono-dimensionale.

Sia A = d2/dx2 con u(x) e v(x) due funzioni definite nell’intervallo [0, L]. L’integrale a primo

membro delle equazioni (32) applicato ad u(x) e v(x) viene sviluppato per parti, ottenendo:

∫ L

0

d2u

dx2v dx =

du

dxv

]L

0

−

∫ L

0

du

dx

dv

dxdx (33)

Affinche il termine a secondo membro della (33) esista e ora sufficiente assicurare che u(x) e v(x)

siano entrambe derivabili una sola volta, mentre la forma iniziale dell’integrale necessitava di una

funzione u(x) derivabile almeno due volte. Quando u(x) = v(x), come nel caso del metodo varia-

zionale di Galerkin, l’applicazione della (33) consente di abbassare il grado di regolarita richiesto

per la funzione base prescelta. In altri termini, le funzioni base possono essere scelte da un sotto-

spazio Fn di ordine inferiore, nel senso che contiene funzioni meno regolari, e quindi piu “semplici”.

Ad esempio, nel caso appena esaminato diventa sufficiente utilizzare funzioni base lineari anziche

paraboliche.

Ogniqualvolta si applichi una integrazione per parti agli integrali variazionali di Galerkin si

ottiene una cosiddetta formulazione debole a cui corrisponde una soluzione debole un. Come prece-

dentemente chiarito, la soluzione e debole nel senso che possiede caratteristiche di regolarita inferiori

a quelle richieste dal problema differenziale di partenza. Si puo, tuttavia, dimostrare che al cresce-

re della dimensione n del sottospazio contenente un la soluzione debole converge ad una soluzione

generalizzata che coincide con la soluzione esatta dell’equazione differenziale data.

L’integrazione per parti puo essere facilmente estesa a problemi bi-dimensionali. Date u(x, y),

v(x, y) e w(x, y), funzioni sufficientemente regolari definite in un dominio piano R, vale la seguente

7

dsdy

dx

n

n

n

x

y

γ γ

γ

γ

12

1

2

RRδ

Figura 2: Orientazione della coordinata curvilinea definita lungo la frontiera ∂R e della normaleuscente.

identita:

∫

R

(

∂u

∂x+∂v

∂y

)

w dx dy = −

∫

R

(

u∂w

∂x+ v∂w

∂y

)

dx dy +

∫

∂R

w (unx + vny) ds (34)

dove s e la coordinata curvilinea definita lungo la frontiera ∂R e nx e ny sono i coseni direttori

della normale n uscente da ∂R (Figura 2). In sostanza, si “scarica” la derivata da u(x, y) e v(x, y)

e la si “carica” su w(x, y) al costo di un integrale aggiuntivo da calcolarsi lungo la frontiera del

dominio R. La formula (34) e nota come lemma di Green e puo essere estesa in modo immediato

ad un dominio tri-dimensionale.

Dimostriamo la relazione (34). Si sviluppi l’integrale delle derivate dei prodotti di u(x, y) e

v(x, y) per w(x, y):

∫

R

[

∂

∂x(uw) +

∂

∂y(vw)

]

dxdy =

∫

R

(

u∂w

∂x+ v∂w

∂y

)

dxdy +

∫

R

(

∂u

∂x+∂v

∂y

)

w dxdy (35)

Il primo membro dell’equazione (35) puo essere sostituito da un integrale calcolato sulla frontiera

∂R usando il teorema di Green in base al quale se P e Q sono funzioni continue e derivabili nella

chiusura di R vale:∫

R

(

∂Q

∂x−∂P

∂y

)

dxdy =

∫

∂R

(P dx+Qdy) (36)

Applicando il risultato (36) alla (35) in cui si e assunto Q = uw e P = −vw si ricava:

∫

∂R

w (u dy − v dx) =

∫

R

(

u∂w

∂x+ v∂w

∂y

)

dxdy +

∫

R

(

∂u

∂x+∂v

∂y

)

w dxdy (37)

Introduciamo ora nell’integrale a primo membro della (37) una coordinata curvilinea s lungo la

frontiera ∂R. Dalla rappresentazione grafica di Figura 2 si osserva che, in base all’orientazione

dell’elemento infinitesimo ds e della normale n, vale:

dy

ds= cos γ2 = nx −

dx

ds= cos γ1 = ny (38)

8

Esplicitando nella (37) l’ultimo integrale ed introducendo le definizioni (38) segue immediatamente

la (34), chiudendo quindi la dimostrazione.

A titolo di esempio, applichiamo il lemma di Green al caso in cui si consideri l’operatore

laplaciano A = ∇2 con due funzioni u(x, y) e v(x, y) definite in R:

∫

R

(

∂2u

∂x2+∂2u

∂y2

)

v dxdy = −

∫

R

(

∂u

∂x

∂v

∂x+∂u

∂y

∂v

∂y

)

dxdy +

∫

∂R

v

(

∂u

∂xnx +

∂u

∂yny

)

ds (39)

E facile osservare che nell’integrale di bordo della (39) compare la derivata normale della funzione

u(x, y):∂u

∂xnx +

∂u

∂yny =

∂u

∂n(40)

Come vedremo piu avanti, questo rende agevole e conveniente il calcolo di tale integrale in quanto

e possibile sfruttare le condizioni al contorno di Neumann. L’utilizzo di una formulazione debole,

quindi, consente di facilitare anche il trattamento delle condizioni al contorno di tipo naturale.

3 Il metodo degli elementi finiti

A seconda della scelta operata per le funzioni base ξ1, ξ2, . . . , ξn si possono ottenere diversi me-

todi variazionali. Il metodo degli elementi finiti e caratterizzato dalla definizione di ξi come un

polinomio di interpolazione a supporto locale continuo a tratti. Questo consente due importanti

vantaggi:

1. i coefficienti αj della combinazione lineare (15) che descrive la soluzione approssimata sono

direttamente i valori della funzione incognita sugli n nodi individuati nel dominio di calcolo

R;

2. la scrittura delle equazioni di Galerkin (32) conduce ad un sistema algebrico di equazioni che

risulta assai sparso, consentendo quindi un importante risparmio di memoria computazionale

con la possibile selezione di un grande numero di nodi su cui calcolare la soluzione numerica

del problema considerato.

3.1 Elementi finiti mono-dimensionali

Sia u(x) una funzione definita nell’intervallo reale [a, b] che risolve un determinato problema dif-

ferenziale del tipo (17)-(18) e siano u1, u2, . . . , un i valori puntuali che essa assume negli n punti

x1, x2, . . . , xn, con x1 = a e xn = b. Un semplice modo per approssimare u(x) in [a, b] consiste nel

collegare ciascuna coppia di punti con un segmento di retta (Figura 3). E intuitivo notare come

al crescere di n la curva spezzata tenda alla funzione esatta, assicurando cosı la convergenza di un

ad u.

9

x x x xa=x2 3

x4

x5

x =bnn−1n−21

Figura 3: Esempio di interpolazione con polinomi lineari continui a tratti.

Il polinomio di primo grado nel j-esimo intervallo fra xj e xj+1 avra la forma:

P (j) (x) = ℓj (x)uj + ℓj+1 (x)uj+1 (41)

dove ℓj e ℓj+1 sono i polinomi base dell’interpolazione di Lagrange a 2 nodi:

ℓj =xj+1 − x

xj+1 − xjℓj+1 =

x− xjxj+1 − xj

(42)

E immediato osservare che ℓj(x) e lineare e vale 1 sul nodo xj e 0 sul nodo xj+1. Analogamente,

ℓj+1(x) e lineare e vale 1 sul nodo xj+1 e 0 sul nodo xj . La funzione un che approssima u(x)

e quindi data dalla somma dei polinomi P (j)(x) definiti in (41) in cui resta inteso che ciascun

polinomio esiste solamente nell’intervallo j-esimo:

un =n∑

j=1

P (j) (x) (43)

Si noti che ogni polinomio P (j)(x) e combinazione lineare dei valori puntuali assunti dalla funzione

u(x) sui nodi selezionati e che polinomi definiti in segmenti adiacenti dipendono entrambi dal

valore assunto da u(x) nel nodo in comune. Raccogliendo tali valori nodali, la funzione un in (43)

puo essere cosı riscritta:

un =n∑

j=1

ξ(1)j uj (44)

dove le ξ(1)j sono risultano funzioni lineari continue a tratti definite come:

ξ(1)1 =

x2 − x

x2 − x1x ∈ [x1, x2]

0 x ∈ ]x2, xn](45)

10

x x x x x xxj−1 j j+1 n−1 n21

ξξ

ξ1

jn

(1)(1)

(1)

1

Figura 4: Comportamento delle funzioni base piramidali.

ξ(1)j =

0 x ∈ [x1, xj−1[x− xj−1xj − xj−1

x ∈ [xj−1, xj ]

xj+1 − x

xj+1 − xjx ∈ ]xj , xj+1]

0 x ∈ ]xj+1, xn]

j = 2, . . . , n− 1 (46)

ξ(1)n =

0 x ∈ [x1, xn−1[x− xn−1xn − xn−1

x ∈ [xn−1, xn](47)

L’equazione (44) risulta formalmente identica alla definizione (15) in cui le funzioni ξ(1)j con j =

1, . . . , n nelle (45)-(47) rappresentano la base dello spazio Fn ed i coefficienti della combinazione

lineare αj coincidono con i valori nodali della funzione incognita. Le ξ(1)j sono anche chiamate

funzioni base piramidali o roof functions. La ragione di tale denominazione appare evidente dal

comportamento riportato in Figura 4. Tali funzioni sono dette a supporto locale perche risultano

non nulle solamente in un intorno del nodo a cui si riferiscono. Si noti che la rappresentazione

di un offerta dalle funzioni base di tipo piramidale e continua ma la sua derivata prima non lo

e sui nodi. Il segmento individuato da una coppia di nodi adiacenti nel quale viene definito un

polinomio interpolatore a supporto locale viene chiamato elemento finito lineare in 1-D.

Qualora la un descritta dalle ξ(1)j non soddisfi ai requisiti minimi di regolarita richiesti dal-

la soluzione di un determinato problema differenziale e possibile migliorare la rappresentazione

aumentando il numero di nodi coinvolti nell’interpolazione locale. Supponiamo, ad esempio, di

considerare delle terne consecutive di nodi. La generica terna individuata dai nodi xj−1, xj e xj+1

e caratterizzata dal seguente polinomio locale:

P (j) (x) = ℓj−1uj−1 + ℓjuj + ℓj+1uj+1 (48)

dove i polinomi base dell’interpolazione di Lagrange a 3 nodi risultano:

ℓj−1 =(xj − x) (xj+1 − x)

(xj − xj−1) (xj+1 − xj−1)ℓj =

(x− xj−1) (xj+1 − x)

(xj − xj−1) (xj+1 − xj)

11

ℓj+1 =(x− xj−1) (x− xj)

(xj+1 − xj−1) (xj+1 − xj)(49)

In sostanza, la funzione u(x) viene approssimata come una successione di archi di parabola. Som-

mando le interpolazioni locali in ciascuna terna di nodi e raccogliendo i valori nodali della funzione

incognita si ottiene la seguente forma per la soluzione approssimata un:

un =

n∑

j=1

ξ(2)j uj (50)

dove le funzioni base ξ(2)j sono polinomi quadratici continui a tratti definiti come:

ξ(2)1 =

(x2 − x) (x3 − x)

(x2 − x1) (x3 − x1)x ∈ [x1, x3]

0 x ∈ ]x3, xn]

(51)

ξ(2)2 =

(x− x1) (x3 − x)

(x2 − x1) (x3 − x2)x ∈ [x1, x3]

0 x ∈ ]x3, xn]

(52)

ξ(2)j−1 =

0 x ∈ [x1, xj−2[

(x− xj−2) (xj − x)

(xj−1 − xj−2) (xj − xj−1)x ∈ [xj−2, xj ]

0 x ∈ ]xj , xn]

j = 4, . . . , n− 2 con passo 2 (53)

ξ(2)j =

0 x ∈ [x1, xj−2[

(x− xj−2) (x− xj−1)

(xj − xj−2) (xj − xj−1)x ∈ [xj−2, xj [

(xj+1 − x) (xj+2 − x)

(xj+1 − xj) (xj+2 − xj−1)x ∈ [xj , xj+2[

0 x ∈ ]xj+2, xn]

j = 4, . . . , n− 2 con passo 2 (54)

ξ(2)n−1 =

0 x ∈ [x1, xn−2[

(x− xn−2) (xn − x)

(xn−1 − xn−2) (xn − xn−1)x ∈ [xn−2, xn]

(55)

ξ(2)n =

0 x ∈ [x1, xn−2[

(x− xn−2) (x− xn−1)

(xn − xn−2) (xn − xn−1)x ∈ [xn−2, xn]

(56)

Una rappresentazione delle funzioni base ξ(2)j nelle (51)-(56) e mostrata in Figura 5. Il segmento

individuato da una terna di nodi consecutivi e chiamato elemento finito quadratico in 1-D.

In modo analogo e possibile aumentare il numero di nodi presenti in ogni elemento finito ed

il corrispondente grado del polinomio base, definendo quindi elementi cubici, quartici, ecc. La

gestione dei polinomi di interpolazione, tuttavia, diventa rapidamente mal-condizionata anche per

un numero relativamente piccolo di nodi per elemento. L’uso di elementi in 1-D con polinomio

interpolatore di grado superiore a 3 e estremamente raro.

12

x x x x x x x xx x x1 2 3 j−2 j−1 j j+1 j+2 n−2 n−1 n

1

ξ

ξ

ξ

ξ

ξ

ξ

1

2

ξj−1

j

j+1

n−1

n

(2)

(2)

(2)

(2)

(2)

(2)

(2)

Figura 5: Comportamento delle funzioni basi quadratiche.

3.2 Elementi finiti lagrangiani

L’interpolazione con polinomi continui a tratti utilizzata per definire gli elementi finiti mono-

dimensionali puo essere facilmente estesa a problemi in 2 o piu dimensioni. Sia u(x, y) la solu-

zione dell’equazione differenziale (17) con condizioni al contorno (18) definita su un dominio R

quadrangolare:

R = (x, y) : a ≤ x ≤ b, c ≤ y ≤ d (57)

Si selezionino n punti x1, x2, . . . , xn nell’intervallo [a, b], con x1 = a e xn = b, ed m punti

y1, y2, . . . , ym nell’intervallo [c, d], con y1 = c e ym = d (Figura 6). Si pensi ora di fissare la

coordinata y, in modo che la funzione u(x, y) dipenda dalla sola variabile x e si possa pertanto

applicare l’interpolazione lineare a tratti (44):

un =n∑

j=1

ξ(1)j (x)u (xj , y) (58)

dove le funzioni base ξ(1)j (x) sono quelle definite in (45)-(47). Consideriamo ora il termine

u(xj , y) nella (58). Quest’ultimo puo a sua volta essere interpretato come una funzione della

sola coordinata y avendo fissato x = xj e pertanto venire anch’esso interpolato mediante la (44):

u (xj , y) ≃m∑

k=1

ξ(1)k (y)u (xj , yk) (59)

dove le funzioni base ξ(1)k (y) sono sempre quelle definite in (45)-(47) in cui la coordinata y

e utilizzata in luogo della x. Introducendo la (59) nella (58) e ponendo per definizione u(xj , yk) =

ujk, si ricava la seguente approssimazione della funzione u(x, y) in R:

unm =

n∑

j=1

m∑

k=1

ξ(1)j (x) ξ

(1)k (y)ujk =

n∑

j=1

m∑

k=1

ξ(1)jk ujk (60)

dove le n×m funzioni base ξ(1)jk sono date dal prodotto delle funzioni base lineari (45)-(47) in x e y.

Per tale ragione esse vengono chiamate anche tensor-product basis functions. Le funzioni base ξ(1)jk

13

x x x x

y

yy

y

y

x1 2 3 n−1 n

1

2

3

m−1

m

a= =b

c=

d=

R

e

j,k j+1,k

j+1,k+1j,k+1

Figura 6: Elemento lagrangiano bilineare a 4 nodi.

sono sempre di tipo polinomiale e a supporto locale, risultando non nulle solamente nella porzione

di piano costituita dai punti (x, y) tali che x ∈ [xj−1, xj+1] e y ∈ [yk−1, yk+1]. Il quadrilatero

definito da due coppie di nodi adiacenti lungo la direzione x e lungo la direzione y (Figura 6)

viene chiamato elemento finito lagrangiano bilineare in quanto le funzioni base sono definite dal

prodotto di polinomi lagrangiani lineari nelle due direzioni coordinate.

Consideriamo un singolo elemento finito bilineare e come in Figura 6. La funzione base ξ(1),ejk

relativa al nodo (xj , yk) e ristretta all’elemento e risulta:

ξ(1),ejk = (a0 + a1x) (b0 + b1y)

= a0b0 + a1b0x+ a0b1y + a1b1xy = c0 + c1x+ c2y + c3xy (61)

cioe un polinomio in x e y il cui termine di grado massimo e quello misto, vale a dire separatamente

lineare in x e y. I coefficienti a0, a1, b0 e b1 sono quelli dei due segmenti di retta che assumono

valore unitario sul nodo (xj , yk) e nullo sui nodi (xj+1, yk) e (xj , yk+1), rispettivamente. I quattro

coefficienti c0, c1, c2 e c3 si possono altresı determinare imponendo quattro condizioni direttamente

al polinomio bilineare ξ(1),ejk . Essendo quest’ultimo una funzione base di tipo interpolante su e le

quattro condizioni sono immediatamente:

ξ(1),ejk (x, y) =

1 se x = xj , y = yk

0 se x = xj , y = yk+1

0 se x = xj+1, y = yk

0 se x = xj+1, y = yk+1

(62)

cioe che ξ(1),ejk valga 1 sul nodo a cui si riferisce e 0 sugli altri. E ovvio che la determinazione di

ξ(1),ejk mediante le (62) oppure passando attraverso il calcolo di a0, a1, b0 e b1 in (61) conduce al

medesimo risultato, vale a dire ξ(1),ejk e unica.

14

x x x x

y

yy

y

y

x1 2 3 n−1 n

1

2

3

m−1

m

a= =b

c=

d=

e

R

j,k+1 j+1,k+1

j+2,k+1

j,k+2 j+1,k+2

j+2,k+2

j+1,k

j+2,kj,k

Figura 7: Elemento lagrangiano biquadratico a 9 nodi.

In modo analogo a quanto fatto per gli elementi finiti mono-dimensionali, e possibile aumentare

il numero di nodi consecutivi presi in ciascuna direzione sui quali eseguire l’interpolazione a tratti.

Ad esempio, utilizzando in ciascuna direzione un’interpolazione quadratica definita su 3 nodi

consecutivi si ricava la seguente rappresentazione per unm:

unm =

n∑

j=1

m∑

k=1

ξ(2)jk ujk (63)

dove ξ(2)jk sono polinomi biquadratici, ovvero separatamente quadratici in x e y, ottenuti mediante

il prodotto delle funzioni base (51)-(56) calcolate lungo x e y. Il quadrangolo definito da 3 terne

consecutive di nodi nelle due direzioni coordinate costituisce un elemento finito lagrangiano bi-

quadratico (Figura 7). Con la notazione adottata nella (61) la funzione base ξ(2),ejk relativa al nodo

(xj , yk) e ristretta all’elemento e risulta (Figura 7):

ξ(2),ejk =

(

a0 + a1x+ a2x2) (

b0 + b1y + b2y2)

= c0 + c1x+ c2y + c3xy + c4x2 + c5y

2 + c6x2y + c7xy

2 + c8x2y2 (64)

dove (a0+ a1x+ a2x2) e (b0+ b1y+ b2y

2) sono i polinomi relativi alle due parabole che assumono

valore unitario sul nodo (xj , yk) e nullo sui nodi (xj+1, yk), (xj+2, yk) e (xj , yk+1), (xj , yk+2),

rispettivamente. Analogamente alle (62), i 9 coefficienti c0, c1, . . . , c8 del polinomio base ξ(2),ejk

in (64) possono equivalentemente determinarsi imponendo che tale polinomio valga 1 sul nodo

(xj , yk) e 0 su tutti gli altri dell’elemento e.

Aumentando il numero di nodi consecutivi utilizzati in ciascuna delle due direzioni coordinate

si possono definire vari elementi finiti lagrangiani le cui funzioni base sono sempre determinate

mediante prodotti delle funzioni base mono-dimensionali. L’estensione degli elementi biquadratici

in Figura 7 porta agli elementi finiti lagrangiani bicubici, costituiti da 16 nodi con funzioni base

15

Figura 8: Elementi serendipity a 8 (sinistra) e 12 (destra) nodi.

separatamente cubiche nelle due variabili x e y. I 16 coefficienti c0, c1, . . . , c15 del polinomio base

per il nodo (xj , yk) si determinano al solito imponendo che ξ(3),ejk valga 1 su tale nodo e 0 su tutti

gli altri. Si potrebbero definire elementi lagrangiani di ordine superiore con un numero di nodi

rapidamente crescente. La determinazione accurata dei polinomi base, tuttavia, puo diventare

facilmente un problema mal-condizionato.

L’estensione degli elementi lagrangiani a domini tri-dimensionali e immediata. Le funzioni base

si ottengono sempre mediante il prodotto di polinomi di interpolazione locale mono-dimensionali.

Gli elementi finiti lagrangiani in 3-D sono comunemente chiamati brick e risultano definiti da 8

nodi per polinomi base separatamente lineari in x, y e z, 27 nodi per polinomi base separatamente

quadratici, 64 nodi per polinomi base separatamente cubici, ecc. Come si nota, in domini tri-

dimensionali il numero di nodi per ciascun elemento cresce molto velocemente con l’ordine del

polinomio di interpolazione.

Riprendiamo ora in considerazione l’elemento lagrangiano biquadratico di Figura 7. Se soppri-

miamo il nodo centrale otteniamo un nuovo elemento a 8 nodi (Figura 8) per il quale il generico

polinomio base ξ(2,s),ejk sara caratterizzato da 8 coefficienti:

ξ(2,s),ejk = c0 + c1x+ c2y + c3xy + c4x

2 + c5y2 + c6x

2y + c7xy2 (65)

determinati imponendo che esso valga 1 sul nodo (xj , yk) e 0 su tutti gli altri nodi dell’elemento.

In sostanza, dal polinomio base (64) dell’elemento biquadratico si elimina il contributo legato al

prodotto x2y2, mantenendo quindi la simmetria di ξ(2,s),ejk rispetto alle due variabili x e y. Il nuovo

elemento finito cosı determinato prende il nome di elemento serendipity a 8 nodi. In maniera

analoga, si possono eliminare i nodi centrali anche in elementi di ordine superiore. Ad esempio, se

da un elemento lagrangiano bicubico si sopprimono i 4 nodi interni si ottiene un nuovo elemento

serendipity a 12 nodi (Figura 8). In questo caso si devono eliminare 4 coefficienti dal polinomio

base corrispondente coincidenti per motivi di simmetria con quelli che moltiplicano i termini x2y2,

x2y3, x3y2 e x3y3.

16

y

x e

i

m

j

Figura 9: Elemento finito triangolare a 3 nodi.

i

j

m

ξi

m

j

1

1

1

ξ

ξ

i

i

j

j

m

m

Figura 10: Funzioni base locali per un elemento finito triangolare a 3 nodi.

3.3 Elementi finiti triangolari

Quando il dominio di integrazione presenta una forma poligonale complessa puo essere utile ricor-

rere ad elementi finiti piu flessibili come quelli di forma triangolare. Si consideri l’elemento finito

a 3 nodi i, j e m in Figura 9. La funzione base ξei relativa al nodo i e ristretta all’elemento e deve

assumere valore unitario su i e nullo su j ed m. Essendo tre sole le condizioni disponibili per la

sua determinazione, ξei avra la seguente forma:

ξei = c0 + c1x+ c2y (66)

vale a dire un polinomio lineare in x e y. Tale espressione descrive un piano nello spazio che

assume valore unitario sul nodo di riferimento e nullo sugli altri due (Figura 10).

Determiniamo i valori dei coefficienti c0, c1 e c2 per la funzione base ξei . Le tre condizioni da

imporre sono esplicitamente:

c0 + c1xi + c2yi = 1

c0 + c1xj + c2yj = 0

c0 + c1xm + c2ym = 0

(67)

17

dove xk e yk sono le coordinate del nodo k, con k = i, j,m. Il sistema lineare di equazioni (67)

puo essere riscritto in modo compatto come:

1 xi yi

1 xj yj

1 xm ym

c0

c1

c2

=

1

0

0

⇒ Cc = e1 (68)

Il vettore soluzione c risulta:

c = C−1e1 (69)

dove la matrice dei coefficienti C e sempre invertibile se i nodi i, j ed m sono non allineati.

Ricordando la formula di Laplace per l’inversa di una matrice regolare si ottiene:

C−1 =1

det (C)

xjym − yjxm xmyi − xiym xiyj − xjyi

yj − ym ym − yi yi − yj

xm − xj xi − xm xj − xi

(70)

Il vettore dei coefficienti c0, c1 e c2 per ξei risulta banalmente pari alla prima colonna di C

−1.

Per determinare ξej e ξem e sufficiente sostituire in (69) al vettore dei termini noti i versori e2 =

(0, 1, 0)T e e3 = (0, 0, 1)T , rispettivamente, per cui i coefficienti di tali funzioni base corrispondono

semplicemente alla seconda e terza colonna di C−1 in (70).

E noto che:

det (C) = 2∆e (71)

dove ∆e rappresenta l’area dell’elemento considerato. La funzione base ξei puo essere quindi

riscritta come:

ξei =ai + bix+ ciy

2∆e(72)

con:

ai = xjym − xmyj bi = yj − ym ci = xm − xj (73)

Analogamente, le funzioni base ξej e ξem avranno la forma (72) con coefficienti:

aj = xmyi − xiym bj = ym − yi cj = xi − xm (74)

am = xiyj − xjyi bm = yi − yj cm = xj − xi (75)

Globalmente, la funzione base ξi associata al nodo i sara data dall’unione delle funzioni base

locali ξei definite su tutti gli elementi triangolari aventi in comune il nodo i (Figura 11). La

conseguente rappresentazione della funzione approssimata un risultera, pertanto, continua ma

non derivabile lungo i lati e sui vertici degli elementi triangolari. Tali funzioni base costituiscono

quindi una naturale estensione delle funzioni base lineari in 1-D (Figura 3).

E possibile aumentare l’ordine dell’approssimazione utilizzata per descrivere un introducendo

un nuovo nodo in corrispondenza al punto medio di ciascun lato di un elemento triangolare. Si

18

x

y

i

ξi

Figura 11: Funzione base relativa ad un nodo con elementi finiti triangolari a 3 nodi.

ottengono cosı degli elementi a 6 nodi nei quali la funzione base ξq,ei relativa al nodo i e ristretta

all’elemento e risulta:

ξq,ei = c0 + c1x+ c2y + c3xy + c4x2 + c5y

2 (76)

cioe un polinomio di grado 2 completo in x e y. E immediato, tuttavia, osservare che, pur utiliz-

zando polinomi di grado superiore al caso di elementi triangolari a 3 nodi, la funzione un basata

sulle ξqj , con j = 1, . . . , n ottenute dalla (76) non possiede derivate prime continue sui lati degli

elementi. Globalmente, quindi, la rappresentazione sara semplicemente continua seppur con una

maggiore accuratezza locale rispetto agli elementi triangolari a 3 nodi.

4 Soluzione agli elementi finiti dell’equazione di Poisson

4.1 Caso 1-D

Si consideri l’equazione di Poisson in 1-D definita nel dominio R costituito dall’intervallo [0, L]:

d2u

dx2= f (x) 0 ≤ x ≤ L (77)

con le seguenti condizioni al contorno:

x (0) = 0(

du

dx

)

x=L

= q(78)

Definiamo una soluzione approssimata un come:

un =

n∑

j=1

αjξj (x) (79)

essendo ξj(x), j = 1, . . . , n, delle opportune funzioni base indipendenti. I coefficienti incogniti αj

vengono calcolati mediante il metodo variazionale di Galerkin, risolvendo cioe il seguente sistema

di equazioni:∫ L

0

d2undx2ξi dx =

∫ L

0

fξi dx i = 1, . . . , n (80)

19

Applichiamo ora alla (80) una formulazione debole grazie alla quale l’integrale a primo membro

diventa:

∫ L

0

d2undx2ξi dx =

dundxξi

]L

0

−

∫ L

0

dundx

dξidxdx

= qξi (L)−

(

dundx

)

x=0

ξi (0)−

∫ L

0

dundx

dξidxdx (81)

Poiche nella (81) compaiono solamente derivate del primo ordine delle funzioni base ξi oppure della

soluzione approssimata un e sufficiente scegliere le ξi come polinomi lineari a tratti. Adottiamo,

pertanto, una suddivisione del dominio in n − 1 elementi finiti lineari mediante la selezione di

n nodi tali che x1 = 0 e xn = L. Supponiamo, per semplicita, che tali nodi siano equidistanti,

cosicche, detta l la dimensione di ciascun elemento finito pari a L/(n− 1), valga:

xi = (i− 1) l i = 1, . . . , n (82)

Le funzioni base ξi assumono la forma analitica (45)-(47), mentre il loro andamento e riportato

in Figura 4. Poiche tali funzioni valgono 1 sul nodo i-esimo e 0 sugli altri, i coefficienti αi della

combinazione lineare (79) coincidono con il valore assunto da un nel punto xi e pertanto verranno

indicati con il simbolo ui. Osserviamo, inoltre, che nella (81) compaiono i valori ξi(0) e ξi(L), i

quali risultano non nulli solamente per i = 1 e i = n, rispettivamente. Il sistema di equazioni (80),

con la formulazione debole (81) e la definizione (79), risulta quindi:

n∑

j=1

(

∫ L

0

dξjdx

dξ1dxdx

)

uj = −

∫ L

0

fξ1 dx−

(

dundx

)

x=0

n∑

j=1

(

∫ L

0

dξjdx

dξidxdx

)

uj = −

∫ L

0

fξi dx i = 2, . . . , n− 1

n∑

j=1

(

∫ L

0

dξjdx

dξndxdx

)

uj = −

∫ L

0

fξn dx+ q

(83)

Definendo:

hij =

∫ L

0

dξjdx

dξidxdx (84)

il sistema di equazioni (83) puo essere riscritto in modo compatto come:

Hu = b (85)

dove H ∈ Rn×n e la matrice dei coefficienti hij , u il vettore delle incognite ui, i = 1, . . . , n e b il

vettore contenente i valori a termine noto nelle (83). La soluzione agli elementi finiti dell’equazione

differenziale (77) si otterra risolvendo il sistema lineare (85).

Valutiamo ora il valore dei coefficienti hij . Consideriamo il caso in cui i = j con i 6= 1 e i 6= n.

Il coefficiente diagonale hii dalla (84) risulta:

hii =

∫ L

0

(

dξidx

)2

dx (86)

20

(a) (b)

i−1 i+1i

ξi

ii−1 i+1i−2 i+2

ξi−1

ξi

ξi+1

(c)

ii−1 i+1 j j+1j−1

ξj

ξi

Figura 12: Assemblaggio dei contributi locali su elementi finiti lineari: (a) termine diagonale;(b) nodi adiacenti; (c) nodi non adiacenti.

in cui la funzione base ξi e quella riportata in Figura 12a. Poiche ξi e non nulla solamente nel-

l’intervallo [xi−1, xi+1], vale a dire nei due elementi finiti che condividono il nodo i, e possibile

restringere il calcolo dell’integrale (86) e scomporlo nella somma di due addendi:

hii =

∫ xi

xi−1

(

dξidx

)2

dx+

∫ xi+1

xi

(

dξidx

)2

dx (87)

In altre parole, il coefficiente hii viene calcolato mediante l’assemblaggio dei contributi locali

valutati sugli elementi finiti che condividono il nodo i. Dalla Figura 12a e dall’espressione analitica

(46) si vede immediatamente che la derivata di ξi vale 1/l nell’intervallo [xi−1, xi] e −1/l in

[xi, xi+1]. Pertanto dalla (87) hii vale:

hii =1

l2

∫ xi

xi−1

dx+1

l2

∫ xi+1

xi

dx

=2

l(88)

Consideriamo ora il caso j = i + 1. Le due funzioni base ξi e ξi+1 sono riportate in Figura

12b. La funzione integranda nella definizione (84) di hij risulta non nulla solo dove ξi e ξj sono

contemporaneamente non nulle, vale a dire nell’elemento finito che collega il nodo i al nodo i+ 1.

L’integrale in (84) puo quindi essere ristretto al solo intervallo [xi, xi+1] dove la derivata di ξi vale

−1/l e la derivata di ξi+1 vale 1/l, ottenendo:

hi,i+1 = −1

l2

∫ xi+1

xi

dx

= −1

l(89)

Anche in questo caso, il coefficiente hi,i+1 e stato ottenuto assemblando i contributi locali cal-

colati su tutti gli elementi finiti che condividono i nodi i ed i + 1, cioe il solo elemento definito

21

dall’intervallo [xi, xi+1]. Con un ragionamento analogo si ricava anche hi,i−1 (Figura 12b):

hi,i−1 = hi,i+1 = −1

l(90)

In Figura 12c e infine riportato il caso in cui j non e un nodo adiacente ad i. Poiche le due

funzioni base ξi e ξj non sono mai contemporaneamente non nulle in nessun intervallo del dominio,

la funzione integranda nella definizione (84) risulta sempre zero e quindi hij = 0. Questo si poteva

altresı dedurre osservando che non esistono elementi finiti che condividono sia il nodo i che il

nodo j. Di conseguenza, la riga i-esima della matrice H possiede solamente 3 coefficienti non

nulli disposti simmetricamente rispetto alla diagonale. La matrice H risulta pertanto simmetrica

e sparsa.

Rimangono da trattare i casi i = 1 ed i = n. Le funzioni base ξ1 e ξn sono non nulle solamente

negli elementi finiti di bordo. L’integrale per il calcolo dei termini diagonali h11 e hnn risulta

pertanto esteso ad un solo elemento finito, dando:

h11 =1

l2

∫ x2

x1

dx =1

l(91)

hnn =1

l2

∫ xn

xn−1

dx =1

l(92)

In definitiva, il sistema di equazioni (85) risulta esplicitamente:

1 −1 0 0 · · · · · · 0

−1 2 −1 0 · · · · · · 0

0 −1 2 −1 · · · · · · 0...

. . ....

0 · · · · · · −1 2 −1 0

0 · · · · · · 0 −1 2 −1

0 · · · · · · 0 0 −1 1

u1

u2

u3...

un−2

un−1

un

=

b1

b2

b3...

bn−2

bn−1

bn

(93)

dove, ricordando il supporto locale delle funzioni base ξi, le componenti del termine noto b sono:

b1 = −l

∫ x2

x1

fξ1 dx− l

(

dundx

)

x=0

bi = −l

∫ xi+1

xi−1

fξi dx i = 2, . . . , n− 1 (94)

bn = −l

∫ xn

xn−1

fξn dx+ lq

La matrice H ricavata in (93), tuttavia, risulta singolare. Cio si verifica immediatamente, ad

esempio, per n = 3. In tal caso, infatti, si ha:

H =

1 −1 0

−1 2 −1

0 −1 1

(95)

22

ove la somma cambiata di segno della prima e terza riga da la seconda riga. Questa situazione

e coerente con il fatto che non sono ancora state imposte in alcun modo le condizioni di Dirichlet

in x = 0 (equazione (78)) e che, quindi, a meno di tali condizioni l’equazione differenziale (77)

ammette un’infinita di soluzioni. Per soddisfare la (78) in x = 0 e necessario imporre:

u1 = 0 (96)

Questo si effettua facilmente modificando la prima equazione del sistema (93) in modo da farla

coincidere con la (96). A tal fine e sufficiente modificare la prima riga della matrice H annullando

il termine extra-diagonale e ponendo b1 = 0:

1 0 0 0 · · · · · · 0

−1 2 −1 0 · · · · · · 0

0 −1 2 −1 · · · · · · 0...

. . ....

0 · · · · · · −1 2 −1 0

0 · · · · · · 0 −1 2 −1

0 · · · · · · 0 0 −1 1

u1

u2

u3...

un−2

un−1

un

=

0

b2

b3...

bn−2

bn−1

bn

(97)

La matrice H e ora non singolare ed il sistema (97) puo essere risolto. Si osservi che l’imposi-

zione della condizione di Dirichlet al nodo 1 semplifica notevolmente il calcolo del termine noto

b1 evitando la stima della derivata della funzione un (equazione (94)). Pertanto a termine noto

sara sufficiente considerare il solo effetto delle condizioni di tipo naturale o di Neumann, trascu-

rando qualsiasi contributo dato lungo la porzione di frontiera sulla quale sono imposte condizioni

principali o di Dirichlet. Si osservi, infine, che in corrispondenza dei nodi di bordo sui quali sia

assegnata una condizione di Neumann nulla il termine noto risultante dalle equazioni (80) non va

modificato. Tale condizione al contorno risulta, pertanto, naturalmente soddisfatta.

4.2 Caso 2-D

Si consideri l’equazione di Poisson in 2-D definita in un generico dominio R sottoinsieme chiuso

di R2:∂2u

∂x2+∂2u

∂y2= f (x, y) (x, y) ∈ R (98)

con condizioni al contorno di tipo principale lungo la porzione di frontiera ∂Ru e di tipo naturale

lungo la porzione ∂Rq:

u (x, y) = u (x, y) ∀ (x, y) ∈ ∂Ru∂u

∂n= q (x, y) ∀ (x, y) ∈ ∂Rq

(99)

23

in cui ∂Ru ∪ ∂Rq = ∂R. Sia un un’approssimazione in R della soluzione alla (98) con condizioni

al contorno (99):

un =

n∑

j=1

αjξj (x, y) (100)

dove ξj , j = 1, . . . , n, sono delle opportune funzioni base. Le equazioni variazionali di Galerkin

applicate al problema considerato offrono:

∫

R

(

∂2un∂x2

+∂2un∂y2

)

ξi dxdy =

∫

R

fξi dxdy i = 1, . . . , n (101)

Utilizzando una formulazione debole per l’integrale a primo membro si ricava:

∫

R

(

∂2un∂x2

+∂2un∂y2

)

ξi dxdy = −

∫

R

(

∂un∂x

∂ξi∂x+∂un∂y

∂ξi∂y

)

dxdy +

∫

∂R

∂un∂nξi ds (102)

Come osservato nel caso del problema di Poisson in 1-D, il contributo calcolato sulla frontiera di R

in (102) puo essere ristretto alla porzione di ∂R su cui sono imposte condizioni di Neumann non

nulle. Introducendo la (102) in (101) e ricordando la definizione (100) di un si ottiene il seguente

sistema di equazioni algebriche nelle incognite αj :

n∑

j=1

[∫

R

(

∂ξi∂x

∂ξj∂x+∂ξi∂y

∂ξj∂y

)

dxdy

]

αj = −

∫

R

fξi dxdy +

∫

∂Rq

qξi ds i = 1, . . . , n (103)

Le equazioni (103) mostrano che le funzioni base ξi devono essere derivabili almeno una volta. E

pertanto possibile utilizzare una discretizzazione con elementi finiti triangolari a 3 nodi le cui fun-

zioni base sono caratterizzate dall’espressione analitica locale (72) e dal comportamento riportato

in Figura 11. Avendo tali funzioni base valore unitario su un nodo della griglia di calcolo e nullo

su tutti gli altri, i coefficienti incogniti αj coincidono con i valori della soluzione approssimata sul

nodo j. Con notazione compatta, il sistema di equazioni (103) si puo esprimere come in (85) in

cui gli elementi hij della matrice H risultano:

hij =

∫

R

(

∂ξi∂x

∂ξj∂x+∂ξi∂y

∂ξj∂y

)

dxdy (104)

Dato il supporto locale delle funzioni base prescelte, costituito dall’unione degli elementi triangolari

che hanno in comune un determinato nodo della griglia di calcolo, il calcolo dell’integrale in (104)

si puo scomporre nella sommatoria dei contributi locali valutati sui singoli triangoli:

hij =∑

e

∫

∆e

(

∂ξ(e)i

∂x

∂ξ(e)j

∂x+∂ξ(e)i

∂y

∂ξ(e)j

∂y

)

dxdy (105)

In particolare, analogamente al caso 1-D, nel calcolo del termine hii la sommatoria in (105)

sara estesa a tutti gli elementi che hanno in comune il nodo i; nel calcolo di hij , con i 6= j,

la sommatoria sara estesa a tutti gli elementi che hanno in comune sia il nodo i che il nodo j.

In altri termini, si tratta dei due elementi che condividono il lato i − j. Se i nodi i e j non si

24

riferiscono a nessun lato della griglia di calcolo il corrispondente valore di hij sara nullo in quanto

non esiste alcuna porzione del dominio R in cui ξi che ξj siano contemporaneamente non nulle.

Di conseguenza, la matrice H sara sparsa in quanto per la riga i-esima gli unici termini non nulli

sono il coefficiente diagonale ed i termini extra-diagonali corrispondenti ai lati che hanno almeno

un estremo in i.

Il contributo locale sull’elemento e al coefficiente hij si puo calcolare assai agevolmente. Ri-

cordando, infatti, la (72) si ottiene:

∫

∆e

(

∂ξ(e)i

∂x

∂ξ(e)j

∂x+∂ξ(e)i

∂y

∂ξ(e)j

∂y

)

dxdy =bibj + cicj4∆e

(106)

dove bi, bj , ci e cj sono funzione delle coordinate dei vertici del triangolo e mediante le (73) e (74).

I contributi a termine noto dipendono, infine, dalla forma analitica assunta dalle funzioni f(x, y)

e q(x, y), per il calcolo dei quali e generalmente conveniente ricorrere a formule di integrazione

numerica di tipo gaussiano. Dopo aver ottenuto la matrice H ed il vettore b e sempre necessario

imporre le condizioni di Dirichlet (99) modificando le equazioni corrispondenti ai nodi che cadono

sul bordo ∂Ru:

ui = u (xi, yi) ∀ (xi, yi) ∈ ∂Ru (107)

Con tali modifiche la matrice H risulta non singolare ed il sistema lineare algebrico corrispondente

puo essere risolto.

25

Esercizi

1. Si consideri l’equazione differenziale:

∂2u

∂x2+∂2u

∂y2= 4− 2x2 − 2y2

nel dominio quadrangolare S = (x, y) : −1 ≤ x ≤ 1,−1 ≤ y ≤ 1 con condizioni al contorno

principali (o di Dirichlet) omogenee. Risolvere l’equazione data mediante il metodo varia-

zionale di Galerkin discretizzando il dominio S mediante 16 elementi lagrangiani bilineari di

lato h = 0.5 (Figura 13) e confrontare la soluzione ottenuta sui nodi della griglia di calcolo

con quella ricavata nella esercitazione n. 3 mediante differenze finite.

y

h

h

x

Figura 13:

2. Dimostrare che la matrice H ottenuta risolvendo un’equazione di Poisson sul dominio S del

precedente esercizio con la discretizzazione ad elementi finiti triangolari rettangoli di Figura

14 coincide con quella ottenuta mediante differenze finite del secondo ordine.

y

h

h

x

1 2 3 4 5

6 8 9 10

11 12 13 14 15

16 18 19 20

21 22 23 24 25

17

7

Figura 14:

26

3. Si consideri l’equazione differenziale:

d2u

dx2−du

dx= f (x)

definita nell’intervallo [0, L] con le condizioni al contorno (78). Determinare la matrice H

ottenuta discretizzando il dominio in (n− 1) elementi finiti lineari con dimensione costante

l = L/(n− 1) ed utilizzando il metodo variazionale di Galerkin. Qual e la differenza con la

matrice H riportata nella (97)?

4. Si consideri l’equazione differenziale:

d2u

dx2+ u+ x = 0

definita dell’intervallo [0, 1] con condizioni principali omogenee sul contorno. Risolvere tale

equazione con il metodo variazionale di Galerkin ed assumendo una soluzione approssimata:

u = α1ξ1 (x) + α2ξ2 (x)

con ξ1 = x(1− x) e ξ2 = x2(1− x). Sapendo che la soluzione esatta e :

u (x) =sin x

sin 1− x

confrontare il valore di u con quello di u sui punti x = 0.25, x = 0.5 e x = 0.75.

27