3. Il metodo di Galerkin-elementi finiti - Elettrotecnica · G. Miano, Appunti del Corso di Modelli...

G. Miano, Appunti del Corso di Modelli Numerici per i Campi, 2009

3. Il metodo di Galerkin-elementi finiti

In questa Nota tratteremo della soluzione numerica di equazioni a derivate parziali scalari con condizioni al contorno attraverso il metodo di Galerkin. Illustreremo poi, come caso particolare, il metodo degli elementi finiti. In particolare, affronteremo il problema della soluzione numerica dell’equazione di Poisson, Helmholtz, diffusione e propagazione con il metodo degli elementi finiti. Tutti questi problemi hanno in comune l’operatore laplaciano, che “governa” il comportamento spaziale delle soluzioni. Per questa ragione introdurremo il metodo di Galerkin ed il metodo degli elementi finiti riferendoci all’equazione di Poisson. Per non appesantire la trattazione considereremo dapprima un problema monodimensionale e poi estenderemo i risultati alle situazioni più generali ed alle altre equazioni. Infine, applicheremo il metodo di Galerkin e degli elementi finiti ad equazioni integrali di Fredholm di tipo scalare. 3.1 Formulazione debole e metodo dei residui pesati

L’equazione !

2u = f P( ) "P #$ (3.1)

è la forma forte dell’equazione di Poisson. Con l’espressione “forma forte” intendiamo il fatto che bisogna cercare una funzione definita nel dominio ! che soddisfi l’equazione differenziale (3.1), quindi deve avere adeguate proprietà di regolarità. La soluzione del problema (3.1) deve essere cercata nello spazio delle funzioni continue, con derivate prime e seconde continue, lo spazio funzionale C 2

![ ] , che è uno spazio di dimensione infinita.

Il metodo delle differenze finite si basa su un’approssimazione discreta dell’operatore laplaciano nei nodi della griglia del dominio di definizione dell’equazione di Poisson e

2


l’equazione di Poisson è imposta nella forma forte. C’è un altro modo di risolvere lo stesso problema, che è allo stesso tempo generale e molto elegante. Invece di approssimare l’operatore laplaciano si approssima lo spazio in cui si cerca la soluzione: si cerca la soluzione dell’equazione di Poisson in uno spazio funzionale a dimensione finita che approssima lo spazio funzionale di definizione del problema. Allo scopo di ridurre il grado di regolarità richiesto alle funzioni di questo spazio il problema sarà riformulato nella cosiddetta forma debole.

Una soluzione approssimata dell’equazione (3.1), con le appropriate condizioni al contorno, è ottenuta cercando la soluzione in un opportuno spazio di funzioni di dimensione finita M , definite su tutto ! fin sulla frontiera, che indichiamo con U

M. Sia

BM

= w1r( ),w2 r( ),...,w

Mr( ){ } (3.2)

una base1 di UM

. Le funzioni di base wk= w

kr( ) possono essere, ad esempio, polinomi,

funzioni trigonometriche, funzioni lineari a tratti, funzioni costanti a tratti. La soluzione approssimata è cercata nella seguente forma

uMr( ) = a

kwkr( )

k=1

M

! . (3.3)

Lo spazio U

M ha la seguente caratteristica: al crescere della dimensione M la funzione

uM

approssima meglio la generica funzione dello spazio funzionale di definizione del problema, ovvero lo spazio funzionale a dimensione finita U

M tende per M !" allo

spazio funzionale di dimensione infinita in cui bisogna cercare la soluzione del problema. Le funzioni w

kr( ) sono dette “funzioni di forma”, mentre i coefficienti a

k sono i

cosiddetti “gradi di libertà”. Allora, la soluzione del problema è ricondotta alla determinazione dei gradi di liberta dell’espansione (3.3). Come si scelgono le funzioni di forma? Come si determinano i gradi di libertà? Ora affronteremo queste due questioni. Affronteremo prima il problema di come determinare i gradi di libertà, poi affronteremo il problema della scelta delle funzioni di forma. 3.1.1 Formulazione debole dell’equazione di Poisson monodimensionale

Come al solito, per introdurre nel modo più semplice possibile l’essenza di un’idea, esemplificheremo dapprima attraverso un problema semplice: determinare la funzione u = u x( ) che verifica l’equazione

d2u

dx2= f x( ) per x ! 0,l( ) , (3.4)

1 Attraverso una opportuna combinazione lineare delle funzioni della base è possibile rappresentare una qualsiasi funzione di U

M.

3


e assegnate condizioni ai limiti, che per il momento non è necessario specificare; u appartiene a C 2

0,l[ ] , ovvero lo spazio delle funzioni continue fino alla deriva seconda definite in 0,l[ ] e f appartiene a C 0

0,l[ ] , ovvero lo spazio delle funzioni continue definite in 0,l[ ] . La (3.4) è la formulazione forte dell’equazione di Poisson monodimensionale. Introduciamo il residuo dell’equazione (3.4) così definito

R w x( ){ } !d2w

dx2" f x( ) per x ! 0,l( ) . (3.5)

Le soluzioni dell’equazione (3.4) sono le funzioni derivabili almeno due volte che annullano il residuo in ogni punto di 0,l( ) ,

R u x( ){ } = 0 per x ! 0,l( ) . (3.6)

Sia r = r x( ) una funzione definita in 0,l[ ] e ivi identicamente nulla,

r x( ) = 0 . (3.7)

Allora, si ha che

p x( )r x( )dx = 00

l

! (3.8)

per ogni funzione p = p x( ) continua definita in 0,l[ ] . E’ vero anche il contrario: se la (3.8) è verificata per ogni funzione continua p = p x( ) definita in 0,l[ ] , allora la funzione r = r x( ) deve essere necessariamente uguale a zero nel dominio 0,l[ ] (ad eccezione di un numero finito di punti). Infatti, se r x( ) fosse diversa da zero in un intervallo di 0,l[ ] di lunghezza diversa da zero esisterebbe almeno una funzione p x( ) per cui la (3.8) non sarebbe verificata. Allora le (3.7) e (3.8) sono equivalenti, cioè

r x( ) = 0 in 0,l[ ]! p x( )r x( )dx = 0 " p x( )0

l

# $C 00,l[ ] . (3.9)

Le funzioni p x( ) sono denominate “funzioni peso”. Questo è il metodo dei residui pesati: imporre la condizione (3.7) equivale a imporre che le “proiezioni di r x( ) lungo ogni elemento dello spazio C 0

0,l[ ]” sia uguale a zero. La (3.7) impone che r x( ) sia uguale a zero in ogni punto del dominio 0,l[ ] , mentre la (3.8) impone che la media di

4


r x( ) pesata per una qualunque funzione peso sia uguale a zero. In realtà questa equivalenza sussiste anche se si restringe lo spazio delle funzioni peso allo spazio delle sole funzioni continue e derivabili definite in 0,l[ ] . Un’ulteriore restrizione che non inficia la validità della (3.9) consiste nel considerare funzioni peso che siano uguali a zero agli estremi dell’intervallo.

Torniamo al nostro problema. Per determinare la soluzione del problema in esame non è necessario imporre che sia verificata la (3.6), ma basta imporre che

p x( )R u x( ){ }dx = 0 ! p x( )0

l

" #C 00,l[ ] . (3.10)

Per non appesantire la notazione è utile introdurre la seguente notazione:

w,v = w x( )v x( )dx0

l

! . (3.11)

w,v rappresenta l’estenzione del prodotto scalare 2 tra due vettori di uno spazio a

dimensione finita alle funzioni continue definite in 0,l[ ] , che costituiscono uno spazio di dimensione infinita. Utilizzando la notazione appena introdotta, la (3.10) può essere così riscritta 3

p, !!u = p, f " p x( )#C 00,l[ ] . (3.12)

Restringendo la (3.12) alle funzioni continue e derivabili C1

0, l[ ] , possiamo applicare la formula dell’integrazione per parti e otteniamo

p, !!u = " !p x( ) !u x( )dx0

l

# + p !u x=

x= l= " !p , !u + p l( ) !u l( ) " p 0( ) !u 0( ) . (3.13)

Allora, determinare le soluzioni dell’equazione (3.4) equivale a determinare le funzioni u tali che

! "p , "u + p l( ) "u l( ) ! p 0( ) "u l( ) = p, f # p x( )$C10,l[ ] . (3.14)

2 Nello spazio funzionale lineare U è definito un prodotto scalare se ad ogni u, v !U è associato un numero reale, denotato usualmente con u, v( ) , tale che: 1. u,u( ) > 0 se u ! 0 e u,u( ) = 0 se u = 0 ; 2. u, v( ) = v,u( ) ; 3. !u, v( ) = ! u, v( ) ; 5. u + v,w( ) = u,w( ) + v,w( ) . Per poter definire il prodotto scalare u, v( ) le funzioni u e v devono essere a quadrato integrabile (non devono essere necesariamente continue).

3 Per allegerire la notazione poniamo !u " du / dx e !!u = d2u / dx

2 .

5


Questa è la forma debole dell’equazione di Poisson monodimensionale. A differenza della (3.4) o della (3.12) nell’equazione (3.14) compaiono solo derivate prime della soluzione, nonché delle funzioni peso: il problema è stato riformulato in modo tale da “indebolire” il vincolo sulle proprietà di regolarità richieste all’incognita del problema. Questa è un’altra ragione per cui si dice che la (3.14) è una “formulazione debole” dell’equazione di Poisson, mentre la (3.4) (o l’equivalente (3.12)) è una “formulazione forte”. Come in seguito vedremo è estremamente conveniente dal punto di vista numerico l’aver indebolito il vincolo sulle proprietà di regolarità richieste alla soluzione. La soluzione della forma debole e le funzioni peso non devono essere necessariamente continue e con derivate prime continue, è sufficiente che siano a quadrato integrabile insieme alla derivate prime 4. L’insieme delle funzioni definite in 0,l[ ] a quadrato integrabile e con derivata prima a quadrato integrabile costituisce lo spazio di Sobolev e lo si indica con H1

0,l[ ] . Osservazione La soluzione della forma forte (3.6) è la soluzione che annulla il residuo in ogni punto del dominio, invece la soluzione della forma debole è la funzione tale da annullare la media del residuo sul dominio di definizione, pesata per qualunque funzione peso. Il risultato fondamentale è che se f x( ) è una funzione continua queste due soluzioni coincidono. La forma forte può essere ottenuta dalla forma debole, basta scegliere come funzioni peso le funzioni impulsive,

p x( ) = ! x " #x( ) con !x " 0,l( ) . (3.15)

♦

3.1.1.1 Problema di Dirichlet

Consideriamo, ora, il problema delle condizioni al contorno. Nel problema di Dirichlet bisogna imporre

u 0( ) = g0 ,

u l( ) = gl .

!"#

(3.16)

4 Le funzioni a quadrato sommabile sul dominio ! sono tutte le funzioni u definite in ! per le quali esiste l’integrale u

2

!" dV ed è limitato. L’insieme di queste funzioni è uno spazio di Hilbert e lo si indica

con L2 !( ) . Le funzioni u con derivata prima a quadrato integrabile sono tutte le funzioni per cui esiste

anche l’integrale !u"#

2

dV ed è limitato. L’insieme di queste funzioni è uno spazio di Sobolev e lo si

indica con H 1!( ) .

6


Data l’arbitrarietà delle funzioni peso, esse possono essere scelte uguali a zero o diverse da zero agli estremi del dominio 0,l[ ] . A seconda di questa scelta cambia la (3.14). E’ evidente, allora, che quando non sono noti i valori della derivata prima agli estremi di 0,l[ ] , come in questo problema, conviene scegliere le funzioni peso in modo tale che

p 0( ) = p l( ) = 0 . (3.17)

In questo caso, la (3.14) diventa

! "p , "u = p, f # p x( )$H0

10,l[ ] (3.18)

dove H

0

10,l[ ] indica il sottospazio di H 1

0,l[ ] che verifica le (3.17). Come in seguito vedremo per altri tipi di problemi conviene fare altri tipi di scelta. La (3.18) deve essere risolta imponendo esplicitamente le condizioni al contorno (3.16). Sia u

!= u

!x( ) una funzione di H 1

0,l[ ] tale che

u!

0( ) = g0, u

!l( ) = gl . (3.19)

La funzione u così definita

u x( ) = u!x( ) + u x( ) (3.20)

verifica le condizioni al contorno (3.16) se la funzione u = u x( ) verifica le condizioni al contorno omogenee

u 0( ) = u l( ) = 0 . (3.21) La funzione u x( ) deve essere tale che

! "p , "u = p, f + "p , "u#

$ p x( )%H0

10,l[ ] . (3.22)

Questa formulazione del problema mette bene in evidenza come le condizioni al contorno contribuiscono al termine noto dell’equazione. 3.1.1.2 Problema di Neumann

Nel problema di Neumann bisogna imporre che

!u 0( ) = q0 ,

!u l( ) = ql .

"#$

(3.23)

7


Sostituendo la (3.23) nella (3.14) otteniamo che

! "p , "u = p, f + p l( )ql ! p 0( )q0 # p x( )$H 1

0,l[ ] . (3.24) Imporre la (3.24) è equivalente ad imporre l’equazione di Poisson in forma debole con le condizioni al contorno (3.23). Infatti, la (3.24) è equivalente a (basta applicare di nuovo l’integrazione per parti)

p, !!u " f( ) = p l( ) !u l( ) " ql#$ %& " p 0( ) !u 0( ) " q0#$ %& ' p x( )(H 10,l[ ] . (3.25)

Siccome la (3.25) deve essere verificata per ogni p x( )!H 1

0,l[ ] , deve essere verificata in particolare per ogni p x( )!H0

10,l[ ] . Imponendo che la (3.25) sia verificata per ogni

p x( )!H0

10,l[ ] si ha che

p, !!u " f( ) = 0 p x( )#H0

10,l[ ] , (3.26)

quindi !!u = f . Utilizzando questo risultato la (3.25) dà

p l( ) !u l( ) " ql#$ %& " p 0( ) !u 0( ) " q0#$ %&=0 ' p x( )(H 10,l[ ] . (3.27)

Imponendo, prima, che la (3.27) sia verificata per p 0( ) ! 0 e p l( ) = 0 , e poi che sia verificata per p 0( ) = 0 e p l( ) ! 0 si ottengono immediatamente le condizioni al contorno (3.23). 3.1.2 Formulazione debole dell’equazione di Poisson in 3D

Il metodo dei residui pesati e la formulazione debole possono essere estesi a problemi nel piano e nello spazio senza richiedere alcuna condizione di rilievo sulla geometria del dominio di definizione. Consideriamo l’equazione

!

2u = f in ! , (3.28)

dove ! è un generico dominio e f è una funzione limitata. Le soluzioni devono appartenere allo spazio C 2

![ ]5, cioè allo spazio delle funzioni derivabili almeno due volte, definite in ! . Introduciamo il residuo

5 Con C p

![ ] indichiamo lo spazio delle funzioni definite in ! derivabili p volte.

8


R u r( ){ } ! "

2u # f . (3.29)

Il problema (3.28) è equivalente a determinare il campo scalare u tale che

p r( )!""" R u r( ){ }dV = 0 #p $C 0 ![ ] . (3.30)

Questa è l’estenzione del metodo dei residui pesati per l’equazione di Poisson tridimensionale: determinare le soluzioni della (3.28) equivale a determinare le funzioni u r( ) per le quali la proiezione del residuo dell’equazione (3.28) lungo ogni elemento dello spazio C 0

![ ] è uguale a zero. Come nel caso monodimensionale, questa equivalenza sussiste anche se si restringe lo spazio delle funzioni peso allo spazio C1

![ ] delle sole funzioni continue e derivabili definite in ! . Un’ulteriore restrizione, che non inficia questa equivalenza, consiste nel considerare funzioni peso uguali a zero sulla frontiera di ! . Consideriamo, ora, la (3.30) per funzioni peso appartenenti a C1

![ ] . Utilizzando l’identità vettoriale

! " p!u( ) = p!2u + !p( ) " !u( ) , (3.31)

e la formula di Green

! " p!u( )dV#$$$ = p !u " n( )dS

%#$$ , (3.32)

dove n è il versore normale alla superficie !" che delimita ! , uscente da ! , la (3.30) diventa

! "p #"udV$%%% + p

&u

&ndS =

&$!%% pfdV

$%%% , (3.33)

dove

!u

!n" n #$u (3.34)

è la derivata normale di u sulla superficie !" . La (3.33) è la “forma debole” dell’equazione di Poisson e la (3.28) è la corrispondente “forma forte”. Posto

a,b ! a r( )"### $b r( )dV (3.35)

9


la (3.30) può essere così riscritta

! "p,"u = p, f ! p#u

#ndS

#$!%% &p 'H 1 $( ) . (3.36)


Consideriamo, ora, il problema delle condizioni al contorno. Nel problema di Dirichlet bisogna imporre che

u

!"= g . (3.37)

Come nel caso monodimensionale, data l’arbitrarietà delle funzioni peso, esse possono essere scelte uguali a zero o diverse da zero su !" . A seconda di questa scelta cambia la (3.36). Quando non sono noti i valori della derivata normale di u su !" conviene scegliere le funzioni peso in modo tale che

p!"

= 0 . (3.38) e la (3.36) diventa

! "p,"u = p, f #p $H0

1%( ) (3.39)

dove H

0

1!( ) indica il sottospazio di H 1

!( ) in cui la (3.38) è verificata. La (3.39) deve essere risolta imponendo esplicitamente le condizioni al contorno (3.37). Sia u

!= u

!r( )

una funzioni di H 1!( ) tale che

u

! !"= g . (3.40)

La funzione u così definita

u r( ) = u!r( ) + u r( ) (3.41)

verifica le condizioni al contorno (3.37) se la funzione u = u r( )

u!"

= 0 . (3.42) La funzione incognita u !H

0

1"( ) deve essere tale che

10


! "p,"u = p, f + "p,"u#

$ p %H0

1&( ) . (3.43)

3.1.2.2 Problema di Neumann

Nel problema di Neumann bisogna imporre che

!u

!n!"

= q (3.44)

Sostituendo la (3.44) nella (3.36) otteniamo che

! "p,"u = p, f ! pqdS#$!%% &p 'H 1 $( ) . (3.45)

Imporre la (3.45) è equivalente ad imporre l’equazione di Poisson in forma debole con le condizioni al contorno (3.44). Analogamente al caso monodimensionale, la (3.45) contiene la condizione di compatibilità tra i dati del problema. Infatti, se si particolarizza la (3.45) per p = 1 si ha

1, f = fdV!""" = qdS

#!!"" . (3.46)

Esercizio 3.7 Si dimostri che imporre la (3.45) è equivalente ad imporre l’equazione di Poisson in forma debole con le condizioni al contorno (3.44).

♦

3.2 Metodo di Galerkin

Un modo semplice e generale per risolvere numericamente uno dei problemi in forma debole descritti nel precedente paragrafo consiste nel:

- cercare la soluzione in un sottospazio a dimensione finita U

M dello spazio di

funzioni a cui la soluzione appartiene; M indica la dimensione dello spazio. - imporre che la forma debole sia verificata per ogni funzione peso appartenente

allo stesso spazio a dimensione finita UM

in cui si cerca la soluzione. Consideriamo la formulazione debole dell’equazione di Poisson (con condizioni al

contorno di tipo Dirichlet o di tipo Neumann)

11


trovare u !U : "p,"u = F p( ) #p !U, (3.47) essendo U un opportuno spazio funzionale e F p( ) il termine noto. Come abbiamo già messo in evidenza non è necessario che U sia lo spazio delle funzioni continue e derivabili, è sufficiente che sia un opportuno spazio di Sobolev. Il metodo di Galerkin per l’approssimazione numerica della (3.47) consiste nel cercare una soluzione approssimata uM!U

M, essendo U

M una famiglia di spazi dipendente da un parametro positivo intero

M , tali che

UM!U, dimU

M= M < " . (3.48)

Il problema approssimato assume allora la forma

trovare uM!U

M: "p

M,"u

M= F p

M( ) #pM !UM

. (3.49)

3.2.1 Problema di Dirichlet La forma debole per l’equazione di Poisson con condizioni al contorno tipo Dirichlet

consiste nel determinare la funzione u = u r( )!H0

1"( ) tale che sia verificata la (3.43).

Sia UM

0 un sottospazio a dimensione finita di H0

1!( ) e

BM

0= w

1r( ),w2 r( ),...,w

Mr( ){ } (3.50)

una base di U

M. Una soluzione approssimata della (3.43) può essere espressa nella

seguente forma

uMr( ) = a

hwhr( )

h=1

M

! . (3.51)

I coefficienti a

h sono determinati imponendo che la (3.43) sia verificata per ogni

p !WM

0 . Le funzioni peso possono essere espresse in forma analoga alla (3.51),

p r( ) = chwh r( )h=1

M

! . (3.52)

Al variare dei coefficienti c

h nel campo reale la (3.52) dà tutte le funzioni appartenenti a

UM

0 . Essendo le funzioni di base linearmente indipendenti, imporre la (3.43) per ogni funzione peso (3.52) equivale ad imporla per ogni funzione base w

h. Sostituendo la

(3.51) nella (3.43) e imponendo che sia verifica per ogni funzione di base wh

otteniamo

12


! "ph ,"wk akk=1

M

# = ph , f + "ph ,"u$ h = 1,2,...,M . (3.53)

Questo è un sistema di M equazioni algebriche lineari nelle M incognite a

1,a

2,...,a

M.

Indichiamo con LM

la matrice M ! M con elementi lhk

, per h,k = 1,2,...,M , così definiti

lh,k

! "wh,"w

k (3.54)

con b

M il vettore di dimensione M dei termini noti con elementi b

h, h = 1,2,...,M , così

definiti

bh = ph , f + !ph ,!u", (3.55)

e con a

M il vettore di dimensione M dei coefficienti incogniti dell’espansione (3.51), a

h,

h = 1,2,...,M . Il sistema (3.53) può essere così riscritto

LMaM= !b

M. (3.56)

La matrice

LM

è simmetrica perché

!wh,!w

k= !w

h"!w

kdV

#$$$ = !w

k"!w

hdV

#$$$ = !w

k,!w

h. (3.57)

Questa proprietà è una diretta conseguenza del fatto che abbiamo scelto le stesse funzioni base per rappresentare la soluzione e le funzioni peso. E’ evidente che se avessimo fatto una diversa scelta non sarebbe verificata la (3.57). Essendo

LM

simmetrica si ha che

aTLMc = c

TLMa (3.58)

per ogni a,c !"M . Una matrice con questa proprietà si dice che è auto-aggiunta. Anche l’operatore laplaciano nello spazio C

0

2!( ) è auto-aggiunto, cioè

s r( )!2

zdV

"### = z r( )!2

sdV

"### (3.59)

per ogni coppia di funzioni s r( ), z r( )!C0

2"( ) . Questa proprietà è una diretta

conseguenza dell’identità (3.31) e del fatto che le funzioni sono uguali a zero sulla frontiera di ! . Lasciamo al Lettore la dimostrazione. Non è per caso che !2 e

LM

siano entrambi auto-aggiunti. Si consideri la forma bilineare

aTLMc . Posto

13


s r( ) = ahwhr( )

h=1

M

! , (3.60)

z r( ) = chwhr( )

h=1

M

! , (3.61)

e tenendo presente che

aTLMc = a

hlhkck

k=1

M

!h=1

M

! = ah"w

h,"w

kck

k=1

M

!h=1

M

! , (3.62)

si ha immediatamente

aTLMc = !s,!z = " s,!

2z . (3.63)

La matrice

LM

è definita positiva, cioè

a

TLMa > 0 !a " 0 . (3.64)

Infatti, dalla (3.63) si ha

aTLMa = !s,!s = !s

2dV

"### . (3.65)

Allora si ha

aTLMa > 0 per !s " 0 e

aTLMa = 0 per !s = 0 , cioè per s r( ) costante in

! . Essendo s r( ) uguale a zero sulla frontiera di ! , le funzioni s r( ) costanti sono solo quelle uguali a zero ovunque in ! . In conclusione,

aTLMa = 0 se e solo se a = 0 ,

altrimenti la forma quadratica aTLMa è positiva.

Lo stesso risultato può essere ottenuto risolvendo il problema di Dirichlet attraverso la formulazione variazionale illustrata in § 1.8.4. Il minimo del funzionale J v( ) è cercato nello spazio a dimensione finita U

M

0 . Si lascia al lettore la dimostrazione. 3.2.2 Problema di Neumann

La forma debole per l’equazione di Poisson con condizioni al contorno tipo Neumann

consiste nel determinare la funzione u = u r( )!H 1"( ) tale che sia verificata la (3.45). Si

procede come nel problema di Dirichlet. La soluzione viene cercata in un sottospazio a dimensione finita U

M di H

1!( ) . Una soluzione approssimata della (3.45) può essere

espressa nella stessa forma della (3.51),

14


uMr( ) = a

hwhr( )

h=1

M

! . (3.66)

I coefficienti a

h sono determinati imponendo che la (3.45) sia verificata per ogni

p !UM

. Le funzioni peso possono essere espresse in forma analoga alla (3.66),

p r( ) = chwh r( )h=1

M

! . (3.67)

Al variare dei coefficienti c

h nel campo reale la (3.67) dà tutte le funzioni appartenenti a

WM

. Come per il problema di Dirichlet, essendo le funzioni di base linearmente indipendenti, imporre la (3.45) per ogni funzione peso (3.67) equivale ad imporla per ogni funzione base w

h. Sostituendo la (3.66) nella (3.45) e imponendo che sia verifica

per ogni funzione di base wh

otteniamo

! "ph ,"wk akk=1

M

# = ph , f ! phqdS$%!&& h = 1,2,...,M , (3.68)

ovvero

LMaM= !b

M (3.69)

dove la matrice

LM

e il vettore incognito a hanno espressioni analoghe a quelle ottenute per il problema di Dirichlet, mentre il vettore dei termini noti b

M ha componenti

bh = ph , f ! phqdS"#!$$ . (3.70)

Anche per questo problema la matrice

LM

è simmetrica, ma non è definita positiva. Nel problema di Neumann accade che la matrice

LM

è semidefinita positiva, cioè

aTLMa ! 0 per a ! 0 . (3.71)

A differenza di quanto accade per il problema di Dirichlet, esiste almeno un vettore a ! 0 tale che aTL

Ma = 0 . Infatti, come per il problema di Dirichlet si ha

aTLMa = !s,!s = !s

2dV

"### . (3.72)

Essendo, in questo caso, s r( ) ! 0 sulla frontiera di ! , la funzione s può essere costante e diversa da zero in ! , quindi esiste almeno un a ! 0 per cui

aTLMa = 0 . Gli elementi

15


di a sono i gradi di libertà dell’espansione (3.66) di una funzione costante in ! . Questo risultato, ovviamente, è valido solo se è possibile rappresentare una funzione costante attraverso la base

BM

. Su questa questione torneremo in seguito quando affronteremo il problema della scelta delle funzioni di base.

Assumiamo, per il momento, che attraverso la base BM

sia possibile rappresentare una funzione costante su ! . Essendo

LM

simmetrica e semidefinita positiva si ha che essa ha almeno un autovalore uguale a zero e, quindi, non è invertibile. E’ un problema? In realtà, come ora vedremo, non c’è nessun problema, perché l’insieme di equazioni (3.68) è un insieme di equazioni linearmente dipendenti e il rango della matrice

LM

è M !1 . Le equazioni indipendenti del sistema (3.68) sono M !1 . Ora dimostriamo questo risultato. Siano gli elementi di a i gradi di libertà della funzione unitaria in ! ,

s r( ) = ah ph r( )h=1

M

! = 1 . (3.73)

Moltiplichiamo ambo i membri dell’equazione (3.68) per a

h e poi sommiamo ambo i

membri di tutte le equazioni così ottenute per h = 1,2,...,M . Abbiamo:

! ah"phh=1

M

# ,"wk akk=1

M

# = ah phh=1

M

# , f ! ah phh=1

M

# qdS$%!&& , (3.74)

ovvero

! "s,"wk akk=1

M

# = s, f ! sqdS$%!&& . (3.75)

A causa della (3.73) e della condizione di compatibilità tra i dati (3.46) la (3.75) è un’identità, quindi le M equazioni (3.68) sono linearmente dipendenti.

Il rango della matrice LM

è M !1 . Si elimini, ad esempio, l’ultima riga della matrice

LM

e si indichi con !LM

la matrice M !1( ) " M !1( ) così ottenuta. Per !LM

si ha

!aT !L

M!a = !!s,!!s = !!s

2dV

"### (3.76)

dove

!s r( ) = !ahwhr( )

h=1

M !1

" . (3.77)

L’aver eliminato un grado di libertà non consente, in generale, di rappresentare più una funzione costante sul dominio ! , a meno che non sia identicamente nulla. Di conseguenza, la matrice

!LM

è definita positiva, quindi è invertibile. In conclusione, il sistema (3.69) ammette infinite soluzioni. La soluzione è unica se si fissa il valore di un grado di libertà nell’espansione (3.66). Il valore di questo grado di libertà non può che

16


essere arbitrario. Ciò è in accordo con il fatto che il problema di Neumann ha un’unica soluzione a meno di una costante additiva arbitraria.

3.2.3 Analisi del metodo di Galerkin Ci soffermeremo, ora, brevemente su queste tre questioni generali:

- esistenza ed unicità della soluzione dell’equazione di Galerkin; - stabilità della soluzione dell’equazione di Galerkin; - convergenza della soluzione dell’equazione di Galerkin alla soluzione esatta.

La soluzione del problema di galerkin esiste ed è unica in conseguenza del fatto che la matrice L

M è definita positiva.

La soluzione del problema di Galerkin è stabile perché l’inversa della matrice LM

è limitata per ogni M (quindi anche per M !" ), allo stesso modo di quanto si verifica nel metodo delle differenze finite. Infine la soluzione del metodo di Galerkin converge alla soluzione del problema esatto se per M !" lo spazio U

M tende allo spazio intero U .

3.3 Metodo degli elementi finiti

La scelta delle funzioni di base è fatta, sostanzialmente, in base a due esigenze contrastanti:

i) le funzioni devono appartenere allo spazio H 1!( ) , perché devono avere

derivate a quadrato integrabili; ii) la complessità e il costo computazionale per calcolare la matrice

LM

e il termine noto b devono essere ridotti al minimo.

Come al solito, solo allo scopo di semplificare la trattazione, tratteremo prima il caso

di domini monodimensionali e poi estenderemo i risultati a domini più generali bidimensionali e tridimensionali.

Il metodo degli elementi finiti consiste:

a) nel partizionare il dominio di definizione del problema ! in domini

“elementari”, detti elementi finiti; b) nello scegliere come funzioni di base funzioni di H 1

!( ) che hanno come supporto solo elementi finiti contigui.

Un “elemento finito” è un dominio di ! con dimensioni molto più piccole della più piccola lunghezza caratteristica su cui la soluzione varia in modo significativo. Attraverso una scelta opportuna della forma e del numero di elementi finiti è possibile partizionare adeguatamente domini di forma complessa. Il “supporto” di una funzione è la parte del

17


suo dominio di definizione dove è diversa da zero (quasi ovunque). Già si intuisce che il metodo degli elementi finiti, a differenza del metodo delle differenze finite, consente di risolvere in maniera adeguata problemi definiti su domini geometricamente complessi. Altre scelte sono ovviamente possibili per le funzioni di base. Ad esempio, i metodi spettrali impiegano come funzioni di base funzioni che hanno come supporto l’intero dominio di definizione e sono ortogonali, cioè

wh,w

k= !

hk. (3.78)

Ad un costo computazionale più elevato se confrontato con quello degli elementi finiti non è detto che corrisponda sempre una maggiore accuratezza. 3.3.1 Domini monodimensionali

L’obiettivo di questa sezione è costruire un sottospazio UM

per domini monodimensionali. Introduciamo una partizione

T

!, in generale non uniforme,

dell’intervallo 0,l[ ] in N +1 sottointervalli Ih= x

h!1, x

h[ ] , h = 1,2,...,N +1 . Il generico sottointervallo I

h è un elemento finito della partizione

T

!, Figura 3.1 . Il pedice ! nel

simbolo T

! sta ad indicare la lunghezza dell’elemento finito più grande. La partizione

T

!

definisce una griglia G! costituita dai nodi x

k, k = 0,1,...,N +1 , e dagli elementi I

h. In

ciascun nodo interno incidono due elementi finiti, mentre in un nodo di frontiera incide un solo elemento finito.

Figura 3.1 Griglia di un dominio monodimensionale.

Il modo più semplice ed efficace di costruire una base per UM

consiste nell’utilizzare funzioni che hanno supporto o in due elementi contigui, interni al dominio di definizione, o in un elemento finito che comprende uno dei due nodi di frontiera. Lo spazio generato da funzioni di base che hanno come supporto un solo elemento finito interno non sarebbe adeguato, perché tutte le funzioni di questo spazio sarebbero uguali a zero nei nodi (ricordiamo che le funzioni devono essere continue). Con questa scelta ad ogni nodo della griglia è associata una funzione di base: la funzione w

h= w

hx( ) associata al nodo h ha

supporto negli elementi finiti che incidono nel nodo h . Siccome la funzione whx( ) deve

essere continua, essa deve essere uguali a zero nei nodi contigui h !1( ) e h +1( ) se h è un nodo interno; se h = 0 w

hx( ) deve essere uguale a zero nel nodo 1; se h = N +1

whx( ) deve essere uguale a zero nel nodo N . Una conseguenza molto interessante di

questa scelta è che gli elementi

18


lhk=

dwh

dx

dwk

dxdx

0

l

! h,k = 1,2,...,N (3.79)

della matrice

LM

corrispondenti a coppie di funzioni di base whx( ) , w

kx( ) con supporti

a intersezione vuota sono uguali a zero, di conseguenza LM

è fortemente sparsa. Il supporto di w

hx( ) è l’unione degli intervalli I

h e I

h+1. Di conseguenza solo i supporti di

wh!1

x( ) e wh+1

x( ) hanno intersezione non vuota con il supporto di whx( ) . Pertanto, nel

caso monodimensionale la matrice LM

è tridiagonale. Il numero di funzioni w

hx( ) è uguale al numero dei nodi interni della griglia più il

numero di nodi di frontiera. Le funzioni whx( ) , h = 0,1,....,N +1 , costituiscono un

insieme di funzioni indipendenti perché esse hanno supporti diversi e, quindi, sono una base per uno sottospazio W

M di dimensione M = N + 2 .

E’ comodo scegliere le funzioni di base whx( ) in modo tale che w

hxh( ) = 1 (elementi

finiti lagrangiani). In conseguenza di ciò si ha

u xh( ) = uh = ah , (3.80)

cioè i gradi di libertà coincidono con i valori della funzione incognita nei nodi della griglia G!

. Allora, per gli elementi finiti lagrangiani si ha, in generale,

whxk( ) = !hk (3.81)

dove !

hk è il simbolo di Kronecher

!hk=

1 se h = k,

0 se h " k.

#$%

(3.82)

Abbiamo fino ad ora definito una serie di proprietà generali della base costruita attraverso gli elementi finiti. Resta ancora da specificare l’andamento delle funzioni di base nei rispettivi supporti. 3.3.1.1 Lo spazio X

!

10,l[ ]

Il modo più semplice ed efficace di costruire funzioni di base appartenenti a H 1

0,l( ) consiste nel considerare funzioni continue e lineari a tratti, vedi Nota 7. Indichiamo con X

!

1 lo spazio composto dalle funzioni continue e lineari a tratti sulla partizione T

!,

Figura 3.1. Le funzioni di base di questo spazio che verificano le condizioni sopra descritte sono illustrate in Figura 3.2 e le espressioni analitiche sono

19


w0 x( ) =x1 ! x

x1 ! x0

per x0 " x " x1 , (3.83)

whx( ) =

x ! xh!1

xh! x

h!1

per xh!1 " x " xh

xh+1 ! x

xh+1 ! xh

per xh" x " x

h+1

0 altrove

#

$

%%%

&

%%%

h = 1,2,...,N , (3.84)

wN +1 x( ) =

x ! xN +1

xN +1 ! xN

per xN" x " x

N +1. (3.85)

Figura 3.2 Funzioni di base dello spazio X

!

10, l[ ] .


Applichiamo quanto appena descritto al problema di Dirichlet. La base BM

0 è data da

BM

0= w

1x( ),w2 x( ),...,w

Nx( ){ } , (3.86)

e la funzione di contorno u

!x( ) è data da

u

!x( ) = u

0w0x( ) + u

lwlx( ) . (3.87)

La dimensione di

BM

0 è N . Gli integrali (3.79) possono essere calcolati analiticamente. In particolare, per una griglia uniforme si ha

20


lhk=

!1

"x per k = h +1 e k = h !1

2

"x per h = k

0 negli altri casi

#

$

%%%

&

%%%

(3.88)

dove !x è la lunghezza dell’elemento finito. Per determinare il vettore dei termini noti bisogna calcolare gli integrali

dwh

dx

dw0

dxdx

0

l

! , dwh

dx

dwl

dxdx

0

l

! , wh x( ) f x( )dx0

l

! h = 1,2,...,N . (3.89)

I primi due integrali sono facilmente valutabili analiticamente. Si ha

dwh

dx

dw0

dxdx

0

l

! = "1

#x$h1

, dw

h

dx

dwl

dxdx

0

l

! = "1

#x$hN

. (3.90)

In generale, il terzo integrale non può essere valutato analiticamente, può essere valutato solo numericamente. Si ha immediatamente che

wh x( ) f x( )dxxh!1

xh+1

" = wh x( ) f x( )dxxh!1

xh

" + wh x( ) f x( )dxxh

xh+1

" . (3.91)

Una valutazione approssimata dei due termini nella (3.91), adeguata al nostro scopo, può essere ottenuta utilizzando il metodo dei trapezi (vedi Nota 7). Applicando il metodo dei trapezi si ha


xh

" = #xwh xh!1( ) f xh!1( ) + wh xh( ) f xh( )

2+O #x2( ) , (3.92)

e

wh x( ) f x( )dxxh

xh+1

! = "xwh xh+1( ) f xh+1( ) + wh xh( ) f xh( )

2+O "x2( ) . (3.93)

Siccome w

hxh!1( ) = wh

xh+1( ) = 0 e w

hxh( ) = 1 , combinando le (3.92) e (3.93), e

trascurando i termini O !x2( ) si ha che

21



xh+1

" # $xfh ; (3.94)

si commette un errore che tende a zero come !x

2 per !x2 " 0 . A questo punto possiamo valutare il termine noto. Utilizzando le (3.20), (3.90) e (3.94), dalla (3.55) si ha

bh = !xfh "u0

!x#h1 "

ul

!x#hN h = 1,2,...,N . (3.95)

Riassumendo abbiamo,

LN=1

!x

2 "1 0 0 ... 0

"1 2 "1 0 ... 0

0 "1 2 "1 ... 0

... ... ... ... ... ...

0 ... 0 "1 2 "1

0 ... 0 0 "1 2

. (3.96)

e

bN = !x f1" u

0/ !x

2, f2,..., fN "1

, fN " ul / !x2T

. (3.97) La matrice

LM

ha inversa limitata. In questo caso è facile verificare che il sistema (3.56) coincide con il sistema trovato con il metodo delle differenze finite. Si può dimostrare che l’errore ! , definito come

! =1

lu" # u

M

2

dx

0

l

$%

&'

(

)*

1/2

, (3.98)

tende a zero come !x2 per !x" 0 . Come al solito con u!

= u!x( ) indichiamo la

soluzione esatta del problema in esame. Invece, l’errore sulla derivata prima della soluzione !" , definito come

!" =1

l!u # $ !u

M

2

dx

0

l

%&

'(

)

*+

1/2

, (3.99)

tende a zero come !x per !x" 0 .

22


3.2.1.2 Problema di Neumann

Consideriamo, ora, il problema di Neumann. La base BM

è data da

BM= w

0x( ),w1 x( ),w2 x( ),...,w

Nx( ),w

N +1x( ){ } . (3.100)

La dimensione di

BM

è N + 2( ) . In questo caso la matrice LM

ha dimensione N + 2( ) . Abbiamo per h,k = 0,1,...,N +1 ,

lhk=

!1

"x per k = h +1 e k = h !1;

1

"x2 ! #

h0 ! #hN +1( ) per h = k;

0 negli altri casi.

$

%

&&&

'

&&&

(3.101)

Il termine noto è dato da

bh =

!x

2f0 " q0 ,

!xfh per 1 # h # N ,

!x

2fN +1 + ql .

$

%

&&

'

&&

(3.102)

Il sistema di equazioni da risolvere è,

1

!x

1 "1 0 0 ... 0

"1 2 "1 0 ... 0

0 "1 2 "1 ... 0

... ... ... ... ... ...

0 ... 0 "1 2 "1

0 ... 0 0 "1 1

a0

a1

a2

...

aN

aN +1

= "

0.5 f0!x " q

0

f1!x

f2!x

...

fN!x

0.5 fN +1!x + ql

. (3.103)

Il sistema (3.103) è diverso da quella che si ottiene con il metodo delle differenze finite. Si noti, inoltre, che la matrice è diversa da quella ottenuta per il problema di Dirichlet. Al secondo ordine in !x le equazioni del sistema (3.103) sono linearmente dipendenti. Infatti, la matrice ha rango N +1( ) e la condizione di compatibilità dei dati dà

0.5!xf0" q

0( ) + fh!xh=1

N

# + 0.5!xfN +1+ ql( ) +O !x

3( ) = 0 . (3.104)

23


Basta eliminare una equazione del sistema (3.103), ad esempio l’ultima, e scegliere ad arbitrio il valore di a

N +1, ad esempio, uguale a zero.

3.3.2 Problemi bidimensionali

Il modo più semplice e generale di estendere il metodo degli elementi finiti (che abbiamo appena applicato alla soluzione del problema di Poisson monodimensionale), a problemi bidimensionali consiste nel partizionare il dominio di definizione del problema attraverso domini elementari a forma di triangolo. Consideriamo un dominio piano ! con frontiera !" a forma di poligonale. Sia

T

! una partizione di una poligonale

attraverso elementi finiti a forma di triangolo, Figura 3.3a; come al solito il parametro ! indica la dimensione caratteristica dell’elemento più grande. Anche se nelle situazioni di interesse la frontiera !" non è, in generale, una poligonale, è sempre possibile rappresentarla adeguatamente con una opportuna triangolazione, Figura 3.3b. La naturale estenzione a domini tridimensionali consiste nel partizionare la regione ! attraverso tetraedri, Figura 3.4.

Figura 3.3 Triangolarizzazione di un dominio a forma di poligonale (a) e di un generico dominio (b).

Indichiamo con P1,P

2,...,P

Ni

i nodi interni della partizione T

! e con

PNi +1,PNi +2 ,...,PNi +N f i nodi di frontiera. Ad ogni nodo della griglia

T

! è associata una

funzione di base, quindi, abbiamo M = Ni + N f funzioni di base w

1,w

2,...,w

M. (3.105)

L’insieme dei nodi interni e di frontiera e dei triangoli definisce la griglia G!.

Il supporto della funzione di base wh

associata al generico nodo Ph è l’unione di tutti i

triangoli di T

! che hanno come vertice il nodo P

h. Per costruzione w

hPh( ) = 1 , mentre

deve essere uguale a zero in tutti gli altri nodi della griglia, compresi quelli dei triangoli che hanno in comune il vertice P

h, perché w

hr( ) deve essere continua, quindi

24


whrPk

( ) = !hk . (3.106) Con la scelta fatta, i gradi di libertà rappresentano i valori dell’incognita nei nodi della griglia,

u rPh

( ) = ah . (3.107)

Figura 3.4 Elemento finito nello spazio tridimensionale a forma di tetraedro.

Come nel caso monodimensionale, la matrice LM

è sparsa. Infatti, lhk

, che ricordiamo è dato da

lhk= !w

h"!w

hdV

#$$$ , (3.108)

è diverso da zero solo se w

h e w

k hanno almeno una parte dei rispettivi supporti in

comune. Ciò si verifica solo se i nodi Ph e P

k sono vertici dello stesso triangolo. La

struttura della matrice dipende da come sono ordinati i nodi. 3.3.2.1 Lo spazio X

!

1"[ ]

Come nel caso monodimensionale, il modo più semplice per generare funzioni

appartenenti a H 1!( ) consiste nell’utilizzare le funzioni continue lineari a tratti sulla

partizione T

!. Indichiamo con X

!

1"[ ] lo spazio composto dalle funzioni continue e

lineari a tratti sulla partizione T

!. Le funzioni di base di questo spazio che verificano le

condizioni sopra descritte sono illustrate in Figura 3.5.

25


Figura 3.5

Figura 3.6

Ora determineremo l’espressione analitica della funzione w

h nel triangolo T

r di cui il

nodo Ph è un vertice, Figura 3.6; indichiamo con P

i e Pj gli altri due vertici del triangolo

Tr. I tre nodi sono ordinati in modo tale che la sequenza P

h, P

i ,Pj sia antioraria. Siano

xh , yh( ) , xi , yi( ) e xj, y

j( ) le coordinate in un sistema di coordinate cartesiane dei tre nodi P

h, P

i e Pj , rispettivamente (vedi Figura 3.6). La restrizione della funzione di base

wh

al triangolo Tr, che indichiamo con w

h

r( ) , ha la seguente espressione: wh

r( )x, y( ) = ah

r( )x + bh

r( )y + ch

r( ) (3.109) dove

ahr( )=1

2Aryi ! yh( ) (3.110)

bhr( )=1

2Arx j ! xh( ) (3.111)

26


chr( )=1

2Arxiyj ! x jyi( ) (3.112)

e

Ar =1

2xiyj ! x jyi( ) + x jyh ! xhyj( ) + xhyi ! xiyh( )"

#$% (3.113)

è l’area dell’elemento. Analoghe espressioni si hanno per w

i

r( ) e wj

r( ) (deve essere

preservata la sequenza antioraria dei tre nodi Ph, P

i e Pj . Per il gradiente di w

h

r( ) si ha

!wh

r( )= xa

h

r( )+ yb

h

r( ) , (3.114) quindi !w

h

r( )"!w

n

r( ) per n = h,i, j è dato da

!wh

r( )"!wn

i

r( )= a

h

r( )an

r( )+ b

h

r( )bn

r( ) (3.115) e

!wh

r( ) "!wn

r( )dS = A

rah

r( )an

r( )+ b

h

r( )bn

r( )( )#$$ . (3.116)

All’elemento l

hh della matrice

LM

contribuiscono tutti i triangoli che hanno come vertice il nodo P

h. Il contributo del generico triangolo T

r è

!w

h

r( ) 2dS = A

rah

r( )2+ b

h

r( )2( )"## . (3.117)

All’elemento l

hk contribuiscono tutti i triangoli che hanno in comune entrambi i vertici

Ph e P

k (se almeno uno dei due nodi è interno ci sono due triangoli che hanno in comune

entrambi i vertici Ph e P

k, altrimenti è un solo triangolo).

Tra gli elementi e i nodi della griglia G! c’è una corrispondenza biunivoca: ad ogni

triangolo sono associati in modo univoco i tre nodi che sono vertici del triangolo; ad ogni nodo sono associati in modo univoco i triangoli che hanno quel nodo come vertice. Questa corrispondenza può essere rappresentata attraverso una matrice di incidenza con un numero di righe uguali al numero di triangoli e tre colonne. La r ! esima riga della matrice è associata al triangolo T

r e i suoi elementi sono i tre numeri interi che

identificano i tre nodi che sono vertici del triangolo. Ad esempio, la matrice di incidenza della griglia riportata in Figura 3.7 è (per convenzione i nodi sono ordinati in senso antiorario)

27


C =

1 2 3

3 2 4

2 5 4

2 6 5

7 6 2

1 7 2

. (3.118)

Figura 3.7 3.3.2.2 Problema di Dirichlet

Applichiamo quanto appena descritto al problema di Dirichlet. La base BM

0 è costituita dalle funzioni di base associate ai nodi interni,

BM

0= w

1x( ),w2 x( ),...,w

Ni

x( ){ } , (3.119) e la funzioni di contorno u

!r( ) è data da

u!r( ) = g Pi( )wi r( )

i=Ni +1

Ni +N f

" . (3.120)

La dimensione della base

BM

0 e della matrice LM

è Ni. La funzione u

!r( ) verifica le

condizioni al contorno nei nodi di frontiera (è come se avessimo interpolato linearmente la funzione g definita sulla frontiera). Il termine noto è dato da

bh = wh r( ) f r( )dS + g Pj( ) !wh "!whdS#$$

j=Ni +1

Ni +N f

%#$$ per h = 1,2,...,N

i. (3.121)

Per calcolare il primo integrale nella (3.121) bisogna ricorrere a metodi numerici analoghi al metodo dei trapezi impiegato per domini monodimensionali. Applicando la proprietà di additività degli integrali, il calcolo del primo integrale è ricondotto alla somma di integrali di superficie definiti su triangoli. In Figura 3.8a è riportata una formula di

28


integrazione del primo ordine nell’area del triangolo, in Figura 3.8b è riportata una formula del secondo del secondo ordine, vedi Nota 7.

Figura 3.8 Due formule di integrazione sul triangolo: formula di ordine 1 (a); formula di ordine 2 (b).

L’errore, come nel caso monodimensionale, può essere stimato in due modi:

! =1

A "( )u# $ u

M

"%%

2

dS&

'(

)

*+

1/2

(3.122)

e

!" =1

A #( )$u% & $u

M

#''

2

dS(

)*

+

,-

1/2

; (3.123)

come al solito con u!

= u!r( ) indichiamo la soluzione esatta del problema in esame e con

A !( ) l’area di ! . L’errore ! tende a zero come !2 per !" 0 , invece !" tende a zero come ! , purché i triangoli della griglia non diventino “indefinitivamente ottusangoli” al crescere del loro numero. 3.3.2.3 Problema di Neumann

Consideriamo, ora, il problema di Neumann. La base BM

è data da

BM = w

1,w

2,...,wNi +N f

{ } . (3.124)

In questo caso la dimensione della base e la dimensione della matrice è M = Ni + N f( ) . Il termine noto è dato da

bh = wh fdS!"" # whqdS

$!!" per h = 1,2,...,N

i. (3.125)

29


Il secondo integrale è un integrale di linea sulla frontiera di ! . Esso può essere determinato numericamente applicando, ad esempio, il metodo dei trapezi. Come nel caso monodimensionale, le equazioni del sistema sono linearmente dipendenti al secondo ordine in ! . Basta eliminare una equazione del sistema, ad esempio l’ultima, e scegliere ad arbitrio il valore di a

M, ad esempio, uguale a zero.

3.4 Equazione di Poisson non lineare Impiegando le tecniche di approssimazione dell’operatore laplaciano che abbiamo introdotto precedentemente è possibile risolvere tutti i problemi richiamati nel § 3.1. In questo paragrafo faremo vedere come si risolve un problema di Poisson non lineare. Solo per esemplificare consideriamo il problema monodimensionale

d

2u

dx2= F u; x( ) in 0,l[ ], (3.126)

u 0( ) = g0 ,

u l( ) = gl .

!"#

(3.127)

Assumiamo che la soluzione esista e sia unica. Per risolvere il problema (3.126), (3.127) con gli elementi finiti bisogna costruire prima la forma debole dell’equazione (3.126) con funzioni peso nulle agli estremi, poi scrivere le corrispondenti equazioni di Galerkin e infine rappresentare le funzioni di base e le funzioni peso nello spazio X

!

10,l[ ] . La forma debole dell’equazione (3.126) è:

!p , !u + p,F u; x( ) = 0 per ogni p !H

0

10,l[ ] . (3.128)

Posto

u x( ) = u!x( ) + a

iwi

i=1

N

" x( ) , (3.129)

dove u

!= u

!x( ) è una funzione appartenente a U

M che verifica le condizioni al

contorno, w1,w

2,...,w

N{ } è una base di UM

0 (vedi § 3.2), a1,a

2,...,a

N{ } sono coefficienti incogniti, e imponendo la (3.128) per ogni p ! w

1,w

2,...,w

N{ } si ottiene l’equazione di Galerkin

LNaN+G

NaN( ) = dN , (3.130)

dove

30


aN= a

1,a

2,...,a

N

T , (3.131) gli elementi della matrice

LN

sono dati dalla (3.54), le componenti della funzione vettoriale non lineare GN

sono date da

GhaN( ) = w

h,F u; x( ) per h = 1,2,...,N (3.132)

e le componenti del vettore d

N sono date da

d

h= !w

h, !u

". (3.133)

Costruiamo lo spazio U

N attraverso elementi finiti, in particolare usiamo lo spazio

X!

10,l[ ] (vedi § 3.3) e utilizziamo elementi lagrangiani. Con questa scelta si ha,

l’espressione della matrice LN

è data dalla (3.96) ,

dh =1

!xg0"h1 + gl"hN( ) per h = 1,2,...,N (3.134)

e

wh,F u; x( ) = w

hx( )F u x( ); x!" #$dx

xh%1

xh+1

& ' (xF ah; x

h( ) per h = 1,2,...,N . (3.135)

L’integrale nella (3.135) è stato calcolato in modo approssimato con il metodo dei trapezi (vedi Nota 7), l’errore tende a zero come !x2 per !x" 0 . Utilizzando le (3.135) per le componenti di GN

otteniamo:

Gha( ) = !xF u

h; x

h( ) per h = 1,2,...,N . (3.136) Come nel caso lineare, il modello discreto ottenuto con le differenze finite coincide con il modello discreto ottenuto con gli elementi finiti al secondo ordine in !x . 3.5 Equazione di Helmholtz In questo paragrafo impiegando le tecniche di approssimazione dell’operatore laplaciano descritte precedentemente affrontiamo il problema della soluzione dell’equazione di Helmholtz. Solo per esemplificare consideriamo un’equazione di Helmholtz con condizioni al contorno di tipo Dirichlet in un dominio piano,

!2u " # u = f r( ) in $, (3.137)

31


u!"

= g . (3.138) Se ! è un autovalore dell’operatore laplaciano la soluzione non è unica. Per risolvere il problema (3.137), (3.138) con gli elementi finiti bisogna costruire prima la forma debole dell’equazione (3.137) con funzioni peso nulle agli estremi, poi scrivere le corrispondenti equazioni di Galerkin e infine rappresentare le funzioni di base e le funzioni peso nello spazio X

!

1"[ ] . La forma debole dell’equazione (3.137) è:

!p,!u + " p,u = # p, f per ogni p !C

0

1"[ ] . (3.139)

Posto

u r( ) = u!r( ) + a

iwi

i=1

N

" r( ) , (3.140)

dove u

!= u

!r( ) è una funzione appartenente a U


contorno, w1,w

2,...,w

N{ } è una base per UM

0 , a1,a

2,...,a

N{ } sono coefficienti incogniti, e imponendo la (3.139) per ogni p ! w

1,w

2,...,w


LNaN+ !C

NaN= "b

N, (3.141)

dove gli elementi di

LN

sono dati dalla (3.54), gli elementi di bN

sono dati dalla (3.55) e gli elementi della matrice

C

N sono dati da

c

h,k= w

h,w

k per h,k = 1,2,...,N . (3.142)

L’equazione (3.141) è un sistema di N equazioni algebriche lineari in altrettante incognite. La matrice

C

N è simmetria e definita positiva. Il sistema (3.141) non ha

un’unica soluzione se !" è un autovalore della matrice C

N

!1LN

. Gli autovalori di questa matrice sono tutti reali e positivi, Appendice 1. Lo spazio U

M

0 può essere costruito attraverso elementi finiti, in particolare attraverso gli elementi dello spazio X

!

1"[ ] , vedi § 3.3. Il modello discreto ottenuto con il metodo

degli elementi finiti ha una struttura diversa da quella del modello ottenuto con il metodo delle differenze finite. 3.6 Equazione di diffusione lineare Ora faremo vedere come impiegando le tecniche di approssimazione dell’operatore laplaciano descritte precedentemente sia possibile risolvere un’equazione di diffusione

32


lineare. Solo allo scopo di esemplificare consideriamo un’equazione di diffusione monodimensionale

!u

!t" D

!2u

!x2= f x;t( ) 0<x < l, t > 0 (3.143)

con condizioni al contorno di tipo Dirichlet

u x = 0;t( ) = g0 t( ),

u x = l;t( ) = gl t( ), (3.144)

e la condizione iniziale

u x;t = 0( ) = u0 x( ) . (3.145) Per la condizione di compatibilità deve essere u

0x = 0( ) = g

0t = 0( ) e

u0x = l( ) = gl t = 0( ) . La soluzione è unica. Per risolvere l’equazione di diffusione

(3.143) con le condizioni al contorno (3.144)con gli elementi finiti bisogna costruire prima la forma debole dell’equazione (3.143) con funzioni peso nulle agli estremi, poi scrivere le corrispondenti equazioni di Galerkin e infine rappresentare le funzioni di base e le funzioni peso nello spazio X

!

10,l[ ] . La forma debole dell’equazione (3.143) è:

p,!u

!t+ D "p , "u = p, f per ogni p !H

0

10,l[ ] . (3.146)

Questa equazione deve essere risolta con la condizione iniziale

p,u x;t = 0( ) = p,u0

per ogni p !H0

10,l[ ] . (3.147)

Posto

u x;t( ) = u!x;t( ) + a

it( )w

i

i=1

N

" x( ) , (3.148)

dove u

!= u

!x;t( ) è una funzione appartenente a U


contorno, w1,w

2,...,w


0 (vedi § 3.2), a1t( ),a2 t( ),...,a

Nt( ){ } sono

funzioni incognite, e imponendo la (3.128) per ogni p ! w1,w

2,...,w


CN

daN

dt+DL

NaN= b

N, (3.149)

33


dove gli elementi della matrice

LN

sono dati dalla (3.54), gli elementi della matrice C

N

sono dati dalla (3.142) e gli elementi del vettore bM

sono dati da

bh t( ) = wh , f ! wh ,"u

"

"t! D #wh , #u

" per h = 1,2,...,N . (3.150)

La (3.149) deve essere risolta con la condizione iniziale

C

Ma t = 0( ) = A (3.151)

dove gli elementi del vettore A sono dati da

A

h= w

h, u

0x( ) ! u" x;t = 0( )#$ %& per h = 1,2,...,N . (3.152)

L’equazione (3.149) è un sistema di N equazioni differenziali lineari del primo ordine in N incognite. Costruiamo lo spazio U

M attraverso elementi finiti, in particolare usiamo lo spazio

X!

10,l[ ] (vedi § 3.3). Con questa scelta si ha a

ht( ) = u x

h;t( ) , l’espressione della matrice

LM

è data dalla (3.96),

chk= !

2

3 se h = k

1

6 se k = h ±1

0 negli altri casi

"

#

$$

%

$$

per h,k = 1,2,...,N , (3.153)

le espressioni delle componenti del vettore b

M sono

bh = !xf xh;t( ) +D

!xg0"!x6

dg0

dt

#$%

&'()h1 +

D

!xg0"!x6

dg0

dt

#$%

&'()hN per h = 1,2,...,N (3.154)

e le espressioni delle componenti del vettore A sono

Ah = !xUh "!x

6g0t = 0( )#h1 + gl t = 0( )#hN$% &' per h = 1,2,...,N (3.155)

e U

h= w

h,u

0x( ) / !x . Il calcolo del primo termine sia dell’espressioni (3.154), che

dell’espressione (3.155) è stato effettuato in modo approssimato con il metodo dei trapezi, l’errore tende a zero come !x2 per !x" 0 (vedi Nota 7).

34


Essendo aht( ) = u x

h;t( ) , un altro modo per imporre la condizione iniziale potrebbe

essere a t = 0( ) = U . (3.156)

E’ evidente che la (3.151) con A dato dalla (3.155) è diversa dalla (3.156). Ciò è conseguenza del fatto che nella (3.156) la condizione iniziale (3.145) è imposta “localmente”, mentre nella (3.151) è imposta “globalmente”. Nel limite !x" 0 le due condizioni tendono alla stessa condizione. 3.7 Equazione di propagazione lineare Ora faremo vedere come impiegando le tecniche di approssimazione dell’operatore laplaciano descritte precedentemente sia possibile risolvere un’equazione di propagazione lineare. Solo allo scopo di esemplificare consideriamo un’equazione di propagazione monodimensionale

!

2u

!t2" c

2 !2u

!x2= f x;t( ) 0<x < l, t > 0 (3.157)


u x = 0;t( ) = g0 t( ),

u x = l;t( ) = gl t( ), (3.158)

e le condizioni iniziali

u x;t = 0( ) = u0 x( ),

!u

!tt=0

= !u0x( ).

(3.159)

Per la condizione di compatibilità deve essere u

0x = 0( ) = g

0t = 0( ) e

u0x = l( ) = gl t = 0( ) . La soluzione è unica.

Esercizio 3.7 Si determini il modello discreto agli elementi finiti dell’equazione di propagazione (3.157) con le condizioni al contorno (3.158) e le condizioni iniziali (3.159).

♦

35


3.8 Equazione di diffusione non lineare Ora faremo vedere come impiegando le tecniche di approssimazione dell’operatore laplaciano descritte precedentemente sia possibile risolvere un’equazione di diffusione non lineare. Solo allo scopo di esemplificare consideriamo un’equazione di diffusione non lineare monodimensionale

!b

!t"!2u

!x2= f x;t( )

b = B u( )

#

$%

&%

0<x < l, t > 0 (3.160)


u x = 0;t( ) = g0 t( ),

u x = l;t( ) = gl t( ), (3.161)

e la condizione iniziale

b x;t = 0( ) = b0 x( ) . (3.162) Per la condizione di compatibilità deve essere

b0x = 0( ) = B g

0t = 0( )!" #$ e

b0x = l( ) = B gl t = 0( )!" #$ . Assumiamo che la soluzione sia unica.

Per risolvere il sistema di equazioni (3.160) con le condizioni al contorno (3.161) con gli elementi finiti bisogna costruire prima la forma debole dell’equazione differenziale del sistema con funzioni peso nulle agli estremi, poi scrivere le corrispondenti equazioni di Galerkin e infine rappresentare le funzioni di base e le funzioni peso nello spazio X

!

10,l[ ] . La forma debole dell’equazione differenziale del sistema (3.160) è:

p,!b

!t+ "p , "u = p, f

p,b = p,B u( )

#

$%

&%

per ogni p !H0

10,l[ ] . (3.163)

Il sistema di equazioni (3.163) deve essere risolto con la condizione iniziale

p,b x;t = 0( ) = p,b0

per ogni p !H0

10,l[ ] . (3.164)

Posto

u x;t( ) = u!x;t( ) + a

it( )w

i

i=1

N

" x( ) , (3.165)

36


b x;t( ) = b! x;t( ) + cit( )w

i

i=1

N

" x( ) (3.166)

dove u

!= u

!x;t( ) e b

!= b

!x;t( ) sono funzioni appartenenti a U

M che verificano le

condizioni al contorno, w1,w

2,...,w


0 (vedi § 3.2), a1t( ),a2 t( ),...,a

Nt( ){ } e c

1t( ),c2 t( ),...,c

Nt( ){ } sono funzioni incognite, imponendo la

(3.128) per ogni p ! w1,w

2,...,w


CN

dcN

dt+ L

NaN= d

N,

CNcN= N a

N( ),

!

"#

$# (3.167)

dove

aNt( ) = a

1t( ),a2 t( ),...,a

Nt( )

T , (3.168)

cNt( ) = c

1t( ),c2 t( ),...,c

Nt( )

T , (3.169) gli elementi della matrice

LN

sono dati dalla (3.54), gli elementi della matrice C

N sono

dati dalla (3.142), le componenti del vettore dN

sono

dh t( ) = wh , f !"b""t

#$%

&'(

! )wh , )u" per h = 1,2,...,N (3.170)

e le componenti della funzione vettoriale

N

N sono

NhaN( ) = w

h,B u! x;t( ) + u

kt( )w

kx( )

k=1

M

"#

$%

&

'( ) w

h,!b!!t

per h = 1,2,...,N . (3.171)

Il sistema (3.167) è un sistema composta da N equazioni differenziali lineari del primo ordine e da N equazioni algebriche non lineari in 2N incognite. Esso deve essere risolto con la condizione iniziale

C

NbNt = 0( ) = B (3.172)

dove gli elementi del vettore B sono dati da

Bh= w

h, b

0x( ) ! b" x;t = 0( )#$ %& per h = 1,2,...,N . (3.173)

37


Il sistema (3.167) come nel modello discreto alle differenze finite può essere ridotto a un sistema di sole equazioni differenziali. Costruiamo lo spazio U


X!

10,l[ ] (vedi § 3.3). Con questa scelta si ha a

ht( ) = u x

h;t( ) , c

ht( ) = b x

h;t( ) ,

l’espressione della matrice LN

è data dalla (3.96), le espressioni degli elementi di C

N

sono date dalle (3.153), le espressioni delle componenti del vettore dN

sono

dh = !xf xh;t( ) +D

!xg0"!x6

dc0

dt

#$%

&'()h1 +

D

!xgl "

!x6

dcN +1

dt

#$%

&'()h1 per h = 1,2,...,N , (3.174)

le espressioni delle componenti della funzione vettoriale

N

NaN( ) sono

NhaM( ) = !xB a

h( ) "!x

6

dc0

dt#h1"!x

6

dcN +1

dt#hN

per h = 1,2,...,N , (3.175)

le espressioni delle componenti del vettore B sono

Bh= !xb

0xh( ) "

!x

6b00( )#

h1"!x

6b0l( )#

hN per h = 1,2,...,N , (3.176)

e

c0t( ) = B g

0t( )!" #$ ,

cN +1

t( ) = B gl t( )!" #$ per h = 1,2,...,N . (3.177) I primi termini nelle espressioni (3.174)-(3.176) sono stati valutati approssimativamente tramite il metodo dei trapezi (vedi Nota 7), l’errore tende a zero come !x2 per !x" 0 . Il problema in esame può essere riformulato in modo tale che il modello discreto agli elementi finiti abbia una sola incognita vettoriale di dimensione N . Ad esempio, sostituendo l’espressione di

u = B

!1b( ) nell’equazione differenziale il sistema (3.167) si

riduce a

!b

!t"

!

!xR b( )

!b

!x

#

$%&

'(= f x;t( ) . (3.178)

dove

R b( ) =dB

!1

db. (3.179)

La forma debole di questa equazione è

p,!b

!t+ "p ,R b( ) "b = p, f p !H

0

10,l[ ] . (3.180)

38


L’equazione di Galerkin corrispondente alla (3.180) è

C

N

dcN

dt+G

NcN( ) = dN (3.181)

dove le espressioni delle componenti del vettore d

M sono

dh t( ) = wh , f !"b""t

#$%

&'(

per h = 1,2,...,M , (3.182)

e le espressioni delle componenti della funzione vettoriale

GNcN( ) sono

GhaM( ) = !w

h,R b( ) !b" x;t( ) + c

kt( ) !w

kx( )

k=1

M

#$

%&

'

() per h = 1,2,...,M . (3.183)

L’equazione (3.181) è un sistema di M equazioni differenziali non lineari del primo ordine in M funzioni incognite. Esso deve essere risolto con la condizione iniziale (3.172). Costruiamo lo spazio U


X!

10,l[ ] (vedi § 3.3). Con questa scelta si ha c

h= b x

h;t( ) , le espressioni degli elementi

della matrice C

N sono date dalle (3.153). le espressioni delle componenti del vettore d

M

sono

dh = !xf x1;t( ) "

1

6

dc0

dt#h1 +

dcN +1

dt#hN

$%&

'()

per h = 1,2,...,M , (3.184)

e le espressioni delle componenti della funzione vettoriale

GNcN( ) sono

Nha( ) = !

ch!1

2"xR c

h!1( ) + R ch( )#$ %& +

ch

2"xR c

h!1( ) + 2R ch( ) + R c

h+1( )#$ %&

!ch+1

2"xR c

h+1( ) + R ch( )#$ %& per h = 1,2,...,N .

(3.185)

Il primo termine nell’espressione (3.184) e l’intera espressione (3.185) sono stati valutati approssimativamente tramite il metodo dei trapezi (vedi Nota 7), l’errore tende a zero come !x2 per !x" 0 .

39


3.9 Equazione di propagazione non lineare Consideriamo l’equazione di propagazione non lineare monodimensionale

!2d

!t 2"!2u

!x2= f x;t( )

d = D u( )

#

$%

&%

0<x < l, t > 0 (3.186)


u x = 0;t( ) = g0 t( ),

u x = l;t( ) = gl t( ), (3.187)

e le condizioni iniziali

d x;t = 0( ) = d0 x( ),

!d

!tt=0

= !d0x( ).

(3.188)

Per la condizione di compatibilità deve essere

d0x = 0( ) = D g

0t = 0( )!" #$ e

d0x = l( ) = B gl t = 0( )!" #$ . Si Assuma che la soluzione sia unica.

Esercizio 3.9 Si determini il modello discreto agli elementi finiti dell’equazione di propagazione non lineare (3.186) con le condizioni al contorno (3.187) e le condizioni iniziali (3.188).

♦ 3.10 Altri metodi di soluzione Oltre al metodo di Galerkin esistono altri metodi di soluzione che si basano sull’imporre che il valore medio del residuo dell’equazione sul dominio di definizione, pesato per una qualsiasi funzione peso, sia uguale a zero. Solo per esemplificare faremo riferimento all’equazione di Poisson con condizioni al contorno di tipo Dirichlet. Posto

uMr( ) = u!

r( ) + ahwhr( )

h=1

M

" in ! , (3.189)

40


dove u! verifica le condizioni al contorno e le funzioni di base w

1r( ),w2 r( ),...,w

Mr( )

sono derivabili almeno due volte e sono uguali a zero sulla frontiera di ! , bisogna cercare i gradi di libertà a

1,a

2,..,a

M in modo tale che il residuo

R uM r( );r{ } = !

2uM " f (3.190)

sia mediamente uguale a zero in ! , pesato per opportune funzioni peso p = p r( ) ,

p r( )R uM r( );r{ }d! = 0!" . (3.191)

Se imponiamo che la (3.191) sia verificata per ogni p ! w

1,w

2,...,w

M{ } otteniamo il metodo di Galerkin. 3.10.1 Metodo dei volumi finiti Il dominio ! è suddiviso in M domini elementari !"

i, i = 1,2,...,M . Consideriamo

le M funzioni peso così definite

pir( ) =

1 se r !"i

0 se r #"i

$%&

per i = 1,2,...,M . (3.192)

Imponendo che la (3.191) sia verificata per ognuna di queste funzioni peso si ottiene un sistema di M equazioni linearmente indipendenti in M incognite. 3.10.2 Metodo della collocazione Si considerino M punti P

i= P

iri( ) i = 1,2,...,M in ! opportunamente distribuiti e

consideriamo le M funzioni peso così definite

pir( ) = ! r " ri( ) per i = 1,2,...,M . (3.193)

Imponendo che la (3.191) sia verificata per ognuna di queste funzioni peso si ottiene il sistema di M equazioni linearmente indipendenti in M incognite

R u r

i( );ri{ } = 0 per i = 1,2,...,M . (3.194) Le (3.194) impongono che il residuo sia uguale a zero negli M punti di collocazione.

41


3.10.3 Metodo dei momenti Nel metodo dei momenti si impone che i momenti del residuo fino all’ordine M !1 siano uguali a zero. Ad esempio, per problemi monodimensionali si ha

pix( ) = x i per i = 0,1,...,M !1 . (3.195)

3.10.4 Metodo dei minimi quadrati In questo metodo i gradi di libertà dell’espansione (3.189) sono scelti in modo tale che il funzionale quadratico

r aM( ) = R

2uMr( );r{ }d!

!" (3.196)

sia minimo. Per la condizione di stazionarietà deve essere

!r

!ah

= R!R

!ah

d" = 0# per i = 1,2,...,M . (3.197)

Il punto stazionario è effettivamente un minimo se la matrice hessiana di r a

M( ) è definita positiva in corrispondenza di esso. 3.11 Equazioni integrali di Fredholm

Il metodo di Galerkin-elementi finiti può essere utilizzato anche per risolvere equazioni integrali. Noi qui considereremo due equazioni integrali scalari descritte nella Nota 1: l’equazione di Fredholm di prima specie

K

1P,Q( )! Q( )dSQ

"#$ = f P( ) !P "#$ , (3.198)

dove !" è una superficie chiusa che delimita la regione di spazio ! , il kernel K

1P,Q( )

è dato da

K1P,Q( ) =

1

4!"0

1

rPQ

, (3.199)

rPQ

è il vettore che parte dal punto sorgente Q e termina nel punto campo P e f = f P( ) è una funzione assegnata su !" ; l’equazione integrale di Fredholm di seconda specie

42


K2P,Q( )! Q( )dSQ

"#$ % ! P( ) = f P( ) !P "#$ , (3.200)

dove il kernel K

2= K

2P,Q( ) è dato da

K2P,Q( ) =

1

2!

n P( ) " rPQ

rPQ

2, (3.201)

rPQ

= rPQ/ r

PQ , n = n P( ) è il versore normale alla superficie !" con la normale

uscente dalla regione ! delimitata da !" , f = !2"0n P( ) #E

0P( ) ed E

0= E

0P( )

assegnato su !" . Il kernel K1 dell’equazione integrale di prima specie (3.198) è

singolare per Q! P mentre il kernel K2 dell’equazione di seconda specie (3.200) è

regolare. Risolviamo l’equazione (3.198) con il metodo di Galerkin. Moltiplichiamo ambo i

membri della (3.198) per una funzione w = w P( ) ed integriamo su tutta la superficie !" ,

dSP

!"# dSQw P( )K1 P,Q( )$ Q( )dSQ!"# = w P( ) f P( )dSP

!"# . (3.202)

Imporre che la (3.202) sia verificata per ogni funzione w !L

2"#( ) equivale ad imporre

che la (3.198) sia verificata su !" . Posto

!NP( ) = a

iui

i=1

N

" P( ) , (3.203)

dove

BN= u

1,u

2,...,u

N{ } è una base pelo spazio a dimensione finita UM

che approssima lo spazio funzionale L2 !"( ) a cui appartiene la soluzione dell’equazione integrale (3.202), a

1,a

2,...,a

N sono coefficienti incogniti, sostituendo la (3.203) nella (3.202) e

imponendo la (3.202) per ogni w ! u1,u

2,...,u

N{ } si ottengono le equazioni di Galerkin per l’equazione integrale (3.198)

aj dSP!"# dSQui P( )K1 P,Q( )uj Q( )dSQ

!"# = ui P( ) f P( )dSP!"#

j=1

N

$ i = 1,2,...,N . (3.204)

Posto a = a

1,a

2,...,a

N

T

!RN , introdotta la matrice L

1( ):R

N! R

N , il cui generico elemento è dato da

43


Lij1( )= dSP

!"# dSQui P( )K1 P,Q( )uj Q( )dSQ!"# , (3.205)

ed il vettore dei termini noti b = b

1,b2,...,b

N

T

!RN il cui generico elemento è dato da

bi = ui P( ) f P( )dSP

!"# , (3.206)

il sistema di equazioni (3.204) può essere riscritto come

L1( )a = b . (3.207)

La matrice L1( ) è simmetrica perché il kernel K

1 è simmetrico rispetto alle variabili P e

Q , cioè K1P,Q( ) = K1 Q,P( ) . Inoltre la matrice L1( ) è definita positiva. Infatti si ha

aTL

1( )a = dSP

!"# $ N P( ) dSQK1 P,Q( )$ N Q( )dSQ!"# (3.208)

dove il secondo membro non è altro che l’energia del campo elettrostatico generato dalla distribuzione di carica !

N.

La base BN

può essere costruita con il metodo degli elementi finiti. La superficie !" è approssimata con una griglia di triangoli. Le funzioni di base possono essere, ad esempio, funzioni costanti a tratti associate al baricentro dei singoli triangoli, funzioni lineari a tratti associate ai vertici dei triangoli, … . La caratteristica fondamentale dell’equazione integrale rispetto all’equazione differenziale è che la matrice che rappresenta l’operatore integrale è piena, invece la matrice che rappresenta l’operatore differenziale è sparsa e bandata. Nel calcolo degli elementi della matrice L

1( ) bisogna prestare particolare attenzione al fatto che il kernel è singolare. È fondamentale calcolare correttamente il contributo della singolarità per avere soluzione approssimate che tendano a zero al crescere del numero degli elementi della griglia con l’ordine dell’elemento finito scelto.

L’equazione (3.200) può essere risolta allo stesso modo. Si ottiene: L

2( )a = b (3.209)

dove

Lij2( )= dSP

!"# dSQui P( )K2P,Q( )uj Q( )dSQ

!"# $ ui P( )uj P( )dSP!"# . (3.210)

La matrice L 2( ) non è simmetrica perché K

2P,Q( ) ! K2

Q,P( ) , inoltre non è nemmeno definita. Utilizzando, come nel primo problema, gli elementi finiti per generare la base approssimata si ottiene una matrice piena. Dalla (3.210) è evidente che la matrice L 2( ) ha

44


la diagonale dominante. Anche la diagonale della matrice L1( ) è dominante grazie al

contributo della singolarità.

3. Il metodo di Galerkin-elementi finiti - Elettrotecnica · G. Miano, Appunti del Corso di Modelli...

Documents

Transcript of 3. Il metodo di Galerkin-elementi finiti - Elettrotecnica · G. Miano, Appunti del Corso di Modelli...