45

CAPITOLO III

METODI NUMERICI: CALCOLO AD

ELEMENTI FINITI

Premessa.

In questo capitolo viene presentata ed analizzata una particolare categoria di

metodi numerici: i Metodi ad Elementi Finiti (F.E.M.). L’Analisi ad Elementi

Finiti (F.E.A.) ha storicamente trovato terreno fertile in campo ingegneristico e

fisico, con riguardo, ad esempio, a problemi di analisi strutturale, conduzione del

calore, conduzione elettrica, campi magnetici, fluidodinamica, plasticità in solidi

non lineari, tutti campi in cui il problema fondamentale è la soluzione di (sistemi

di) equazioni differenziali di varia natura. In linea di principio si tratta di

problemi risolvibili con Metodi alle Differenze Finite (F.D.M.) e in effetti, per

certi aspetti, F.E.M. e F.D.M. possono essere considerati “parenti”. Cosa ha

decretato dunque il successo dei F.E.M., considerata anche la loro maggior

complicatezza rispetto ai F.D.M.? Essenzialmente la loro elevata adattabilità a

problemi particolarmente complessi e la capacità di fornire risultati notevolmente

precisi. Punto cruciale nella soluzione di un problema differenziale è il modo in

cui vengono trattate le condizioni al bordo e temporali, e il compito diventa tanto

46

più difficile quanto più complesse sono queste condizioni. Le F.D. si adattano

bene a problemi relativamente semplici e con condizioni regolari, ma mostrano

presto i loro limiti all’aumentare della complessità del sistema analizzato, e ciò

tendenzialmente accade quanto più si cerca di avvicinare la modellizzazione alla

realtà. Non è un caso, quindi, se uno dei primi e più consistenti impieghi dei

F.E.M. fu in ambito aeronautico per la progettazione di profili alari (Clough per

Boeing, 1960). I F.E.M., come i F.D.M., si basano su metodi di discretizzazione,

tuttavia i primi hanno sui secondi l’ulteriore vantaggio di poter utilizzare passi di

discretizzazione variabili. In questo modo l’efficienza di calcolo aumenta

notevolmente, poiché è possibile infittire la mesh solo nei punti di maggior

interesse, lasciando invece maglie più larghe nel resto del domino.

L’interpolazione tra due nodi della mesh può inoltre essere eseguita con varie

tipologie di funzioni, e non esclusivamente lineari, come nel caso delle differenze

finite. Tuttavia anche i F.E.M. hanno, come tutte le rose, le loro spine: la

complessità del metodo è decisamente elevata e questo si riflette completamente

nella difficoltà di programmazione.

Dal momento che il prezzo di un’opzione è descritto da un’equazione

differenziale, il problema della sua individuazione rientra a pieno diritto tra quelli

a cui i F.E.M. possono dare soluzione. In effetti l’implementazione ad elementi

finiti di un metodo di calcolo per il prezzo di un’opzione può sembrare fatica

sprecata, ma in realtà non è così (e comunque non lo sarebbe più di un metodo a

differenze finite). Le opzioni americane, ad esempio, sono descritte da equazioni

differenziali con condizioni al contorno libere e le asiatiche hanno condizioni al

contorno particolarmente complesse: in questi casi i F.E.M. sono di notevole

utilità. La loro applicazione ad opzioni europee è invece giustificata dalla

necessità di testare la validità del metodo: il confronto con la soluzione esatta di

Black-Scholes consente una valutazione diretta e precisa.

47

3.1 Gli Elementi Finiti: generalità.

3.1.1 La configurazione della mesh.

Nei F.D.M. il nucleo centrale dell’analisi è la riformulazione delle derivate

presenti nell’equazione di riferimento tramite espansione in serie, mentre

attenzione praticamente nulla è posta sulla discretizzazione del dominio del

problema: tutto si riduce alla scelta degli intervalli xδ e δτ . Questi possono poi

essere cambiati in funzione della precisione che si desidera raggiungere, ma non

influenzano in alcun modo il metodo di calcolo. Ciò non vale per gli elementi

finiti, in cui la configurazione dell’elemento (cioè l’unità di base in cui è

scomposto il dominio) è fase essenziale. Due sono i motivi principali:

• l’ elemento non ha solo i nodi esterni (cioè gli estremi), ma pu ò averne di

interni, determinando conseguentemente il tipo di funzione interpolante

utilizzato (lineare, quadratica, cubica, …);

• se si intende analizzare con maggior accuratezza un determinato punto del

domino, è necessario infittire la mesh in sua corrispondenza, utilizzando

quindi elementi più piccoli.

Fattori di grande importanza, e che verranno adeguatamente trattati nel paragrafo

relativo alla configurazione degli elementi, sono la condizione di completezza e

la condizione di continuità: qui è appena il caso di ricordare che riguardano le

funzioni di interpolazione e che sono fondamentali per una corretta soluzione del

problema.

ax bx ax cx dx bx

Figura 3.1: esempio di elemento con soli nodi esterni (a) e di elemento con nodi interni (b).

(a) (b)

48

3.1.2 Le tipologie di F.E.M..

Il passo successivo alla definizione della mesh è la scelta della procedura di

calcolo. Diversamente dal caso delle differenze finite, non si opera direttamente

sulla funzione di riferimento, ma su una sua formulazione integrale che può

essere ottenuta tramite due differenti procedure: residui pesati o formulazione

variazionale. Ciò che si ottiene, è, in entrambi i casi, una soluzione approssimata

del problema originale.

I metodi a residui pesati partono dalla considerazione della (generica) equazione

originale ( )( ) QxL =ϕ , dove L è un operatore differenziale agente su ( )xϕ , la

variabile x è definita sul dominio Ω , e Q è il termine noto; si suppone poi che

esista una generica soluzione approssimata ( ) ( ) ( )xaxaax NNψψϕ ++= ...;~11 1 del

problema e la si sostituisce nell’equazione. Dal momento che la soluzione è solo

approssimata, sarà in grado di soddisfare solo in parte la relazione iniziale e si

potrà allora scrivere: ( ) ( ) 0;~ ≠=− axRQL ϕ , dove il termine R indica il residuo (o

errore) derivante dalla sostituzione. Non è possibile far sì che il termine d’errore

scompaia in tutti i punti del dominio (altrimenti ϕ~ sarebbe la soluzione esatta de

problema), ma è possibile richiedere che la somma pesata di tutti gli errori sia

nulla:

( ) ( ) ( ) 0;; == ∫Ω

dx xW axRaxS , (3.1)

1I coefficienti ia sono detti gradi di libertà (D.O.F.) del sistema.

La forma funzionale delle iψ può essere scelta liberamente, tuttavia generalmente si preferisce

una forma polinomiale, che meglio si presta al tipo di calcoli più frequentemente effettuati in

un’analisi ad elementi finiti: derivate ed integrali. Si può considerare

( ) ....;~ 12321

−++++= NN xaxaxaaxaϕ

49

dove ( )xW rappresenta differenti funzioni-peso. La funzione (3.1) è alla base di

tutti i metodi dei residui pesati, che differiscono tra loro solo per il tipo di

funzione-peso utilizzata.

In base al tipo di funzione-peso si distinguono i seguenti metodi:

1. Collocazione.

Per ogni parametro incognito ia si sceglie un punto ix (detto di

collocazione) del domino, e in sua corrispondenza si impone che il residuo

sia nullo:

( )

( ) 0;

0; 1

=

=

axR

axR

N

Si giunge così ad un sistema di N equazioni in N incognite da cui si

ricavano i parametri ia per la soluzione approssimata ( )ax;~ϕ .

Ricordando la misura di Dirac introdotta nel capitolo precedente, può

essere riscritto come:

( ) ( ) ( ) 0; ; ==−∫ axRdxxxaxR i

x

xi

b

a

δ ,

con Ni ,...,1= e bia xxx ≤≤ .

La funzione-peso è data da

( ) ( )ii xxxW −= δ .

50

2. Sottodominio.

Per ogni parametro incognito ia si sceglie un intervallo ix∆ all’interno del

dominio e si impone che la media del residuo in ogni intervallo sia nulla:

( )

( ) 0 ;1

0 ;1

11

=∆

=∆

∫

∫

∆

∆

NxN

x

dxaxRx

dxaxRx

Si ottengono così i parametri ia .

La funzione-peso è

( ) =xWi ( )ii xx −+1 .2

3. Minimi quadrati.

Si minimizza rispetto a ciascun ia l’integrale del quadrato del residuo:

( )

( ) 0 ;

0 ;

2

2

1

=∂

∂

=∂∂

∫

∫

Ω

Ω

dxaxRa

dxaxRa

N

, cioè

( ) ( )

( ) ( )0

;;

0;

;1

=∂

∂

=∂

∂

∫

∫

Ω

Ω

dxa

axRaxR

dxa

axRaxR

N

La funzione-peso è

( ) ( )i

i aaxR

xW∂

∂= ; .

2 E’ la funzione a gradino, di altezza 1 e ampiezza x∆ .

51

4. Galerkin.

Richiede che per ogni parametro ia sia nulla la media dei residui, pesati

tramite le funzioni ( ) ii xW ψ= :

( )

( )∫

∫

Ω

Ω

=

=

0 ;

0 ; 1

dxaxR

dxaxR

Nψ

ψ

La formulazione variazionale affonda le proprie radici nel calcolo delle

variazioni, sviluppato attorno alla metà del 1700 da Eulero. Essa si basa sulla

trasformazione della funzione originale in una sua forma equivalente, ossia la cui

soluzione esatta coincide con quella esatta della formulazione originaria.

L’applicazione di questo metodo è fortemente limitata dall’impossibilità di

trovare una forma variazionale per molti problemi non lineari, mentre è sempre

possibile effettuarne una formulazione a residui pesati. In aggiunta, la soluzione

di una problema variazionale è equivalente a quella ottenibile con i residui pesati

con il metodo Galerkin. La formulazione variazionale è la seguente:

( ) 0=ϕδI ,

dove

( ) ( )( )∫Ω Ω

−+−=

dxd

xdx QLI ϕϕϕϕϕ . (3.2)

L’operatore δ , agendo sul funzionale ( )ϕI , è fonte di piccole variazioni della

funzione ( )xϕ in un intorno della soluzione esatta. Sostituendo nella (3.2) la

52

soluzione approssimata ϕ~ , può essere calcolato l’integrale rispetto a x , poiché le

( )xiψ hanno una forma specifica, e I diventa una normale funzione di parametri:

( )( ) ( )aIaxI =;~ϕ .

( )aI può essere resa stazionaria applicando la condizione 0=dI , cioè

01

=∂∂= ∑

=

N

ii

i

daaI

dI .

Poiché i coefficienti ia possono variare indipendentemente l’uno dall’altro, i

termini della sommatoria devono annullarsi separatamente: si ricavano così i

parametri ia .

Per la soluzione del problema di prezzaggio di un’opzio ne, è stato scelto il

metodo di Galerkin. Innanzitutto i metodi a residui pesati hanno una

connotazione maggiormente applicativa (quindi offrono maggior efficienza in

fase di implementazione), mentre il metodo variazionale si adatta maggiormente

ad approfondimenti teorici (ma non si presta altrettanto bene a una trasposizione

software). In secondo luogo, il metodo di Galerkin è il più efficiente tra quelli a

residui pesati e fornisce comunque risultati equivalenti alla formulazione

variazionale.

3.2 Configurazione degli elementi di classe 0C .

In questo paragrafo si mostra la corretta costruzione di un elemento (inizialmente

considerandone il caso più semplice, e poi casi più complessi e generali) e il

metodo d’assemblaggio tra i vari elementi.

53

Si consideri il dominio Ω e lo si partizioni in N elementi della medesima

lunghezza (che, al momento, non è rilevante): si otterranno 1+N nodi ix ,

11 +≤≤ Ni :

1x 2x 3x 1−Nx Nx 1+Nx

E’ necessario scegliere una particolare forma per la generica soluzione

approssimante ( )ax;~ϕ , nella fattispecie e per le ragioni spiegate nel precedente

paragrafo, la forma ( ) xx 21;~ αααϕ += . Dal momento che questa deve essere la

funzione interpolante per il generico elemento delimitato dai nodi ax e bx , è

necessario assicurare che i coefficienti 1α e 2α siano in corrispondenza dei due

estremi.

0 ax bx

ax bx

elemento 1 elemento N

1 Figura 3.2: Il disegno (a) mostra una forma inappropriata per la costruzione della generica soluzione (in forma di serie di potenze), in cui i coefficienti 1α e 2α non sono in corrispondenza degli estremi dell’elemento da interpolare. Il disegno (b) mostra la corretta funzione interpolante (polinomiale) riferita al medesimo elemento.

( )αϕ ;~ x

( )ax;~ϕ

1a 2a

1α

2α

disegno (a)

disegno (b)

54

Per far ciò, si normalizzano 1α e 2α rispetto a due coefficienti 1a e 2a , scelti in

modo che

( )( ) 2

1

;~;~

aax

aax

b

a

==

ϕϕ

.

Volendo riscrivere la ( )ax;~ϕ in funzione dei nuovi coefficienti, si risolve il

sistema

=+=+

221

121

ax

ax

b

a

αααα

,

da cui risulta

−−

=

−−

=

ab

ab

ab

xxaa

xxaxax

122

211

α

α.

La ( )ax;~ϕ può allora essere riscritta come

( ) xxxaa

xxaxax

axabab

ab

−−+

−−

= 1221;~ϕ ,

o in forma più generale, dopo aver raccolto i coefficienti ia , come

( ) ( ) ( ) ( )∑=

=+=2

12211;~

iii xaxaxaax ψψψϕ , (3.4)

dove

55

( )

( )

−−

=

−−

=

ab

a

ab

b

xxxx

x

xxxx

x

2

1

ψ

ψ. (3.5)

Si noti la seguente proprietà di cui godono le funzioni appena costruite:

( )( )

( )( )

==

==

1

0

0

1

2

2

1

1

b

a

b

a

x

x

x

x

ψψ

ψψ

.

La (3.4) rappresenta la parte di soluzione afferente al generico elemento di

estremi ax e bx (si veda anche la figura 3.3).

ax bx

L’elemento -tipo sinora considerato è caratterizzato dalla presenza di due soli

nodi, posizionati ai suoi estremi: ciò permette di utilizzare solamente funzioni di

interpolazione lineari. Volendo migliorare l’accuratezza della soluzione, due

sono le strade percorribili: ridurre le dimensioni degli elementi e incrementarne il

numero oppure aumentare il grado dei polinomi approssimanti, in modo da

1a

2a

1 1

22ψa

11ψa

1ψ

2ψ

( )ax;~ϕ

Figura 3.3: La soluzione ( )ax;~ϕ come somma delle funzioni-base 1ψ

e 2ψ , pesate secondo i corrispondenti coefficienti 1a e 2a .

56

ottenere curve più complesse. La prima soluzione permette di mantenere i calcoli

a un livello di difficoltà decisamente basso, ma è poco efficiente; la seconda, al

contrario, è molto più efficiente e permette di ottenere la medesima accuratezza

di risultato con un numero di elementi molto inferiore, ma al prezzo di una

struttura numerica più complessa. Nella pratica, la soluzione adottata più

frequentemente prevede un aumento sia del numero degli elementi, sia del grado

del polinomio interpolante, in modo che gli svantaggi dei due metodi si bilancino

reciprocamente.

Per incrementare in grado del polinomio, è necessario aumentare il numero dei

nodi dell’elemento: i nuovi nodi vengono aggiunti internamente agli estremi ed

equispaziati tra loro.

La costruzione dell’elemento procede esattamente secondo la medesima logica

sottostante al caso di interpolante lineare, per cui sembra lecito attendersi, come

risultato finale, una generalizzazione della (3.4):

( ) ( )∑=

=n

iii xaax

1

;~ ψϕ , (3.6)

dove n indica il numero di nodi relativi al singolo elemento. Rispetto al caso

lineare cambiano però le funzioni ( )xiψ , che ora devono tener conto dei nuovi

nodi presenti. Mantenendo la struttura polinomiale per la ( )αϕ ;~ x , questa diventa

( ) ∑=

−− =+++=n

i

ii

nn xxxx

1

1121 ...;~ αααααϕ . (3.7)

Anche in questo caso è necessario normalizzare gli iα rispetto a dei convenienti

parametri ia : per farlo si risolve rispetto agli iα il sistema

57

=

=

∑

∑

=

−

=

−

n

in

ii

n

i

ii

ax

ax

1

1

11

1

α

α

.

Dopo aver inserito nella (3.7) le relazioni appena trovate per gli iα , si effettuano

i raccoglimenti a fattor comune rispetto agli ia , ottenendo così propri

l’espressione estesa della “sospetta” (3.6). La funzione iψ relativa al generico

nodo k è ottenuta mediante il medesimo procedimento del caso

monodimensionale ed è così definita:

( )( )

( )∏

∏

≠=

≠=

−

−

=n

kii

ik

n

kii

i

k

xx

xx

x

1

1

ψ . (3.8)

Le funzioni iψ sono di grado 1−n : ad esempio, se il numero di nodi di un

elemento è 3=n , esso sarà interpolato da tre funzioni quadratiche, ognuna pesata

dal proprio ia . Si noti che per 2=n si ricade, come è lecito aspettarsi, nel caso di

interpolazione lineare analizzato all’inizio del paragrafo.

58

3.3 Costruzione delle matrici locali degli elementi e assemblaggio della

matrice globale.

3.3.1 Costruzione delle matrici locali.

Nel paragrafo (3.1) era stato presentato il metodo dei residui pesati di Galerkin

per un generico problema differenziale ( )( ) QxL =ϕ : ora verrà applicato

1x

1x

1x

1x

2x

2x

2x

2x

3x

3x

3x

3x

1

1

1

1a 3a

( )x3ψ

( )x2ψ

( )x1ψ

( )ax;~ϕ

Figura 3.4: Le tre funzioni-base per una funzione di interpolazione 0C -quadratica. Dall’alto verso il basso si riconoscono le funzioni

( )( )

( )( )3121

321 xxxx

xxxx−−

−−=ψ ,

( )( )

( )( )3212

312 xxxx

xxxx−−

−−=ψ ,

( )( )

( )( )2313

213 xxxx

xxxx−−

−−=ψ .

Nell’ultimo disegno è rappresentata la funzione di interpolazione, data dalla somma pesata delle funzioni-base:

( ) ( ) ( ) ( )xaxaxaax 332211;~ ψψψϕ ++= . Questi sono gli elementi che verranno utilizzati per la soluzione del modello di Black-Scholes.

59

all’equazione di Black -Scholes trasformata 3, considerata su un dominio finito di

x .

( ) ( ) ( )( )( ) 0,lim

0,max0,..

0,,0

11

2

2

21

21

=−=

>+<<−=∂∂−

∂∂

±→

+−

τ

ττ

xu

ee xu ts

LxL xuu

Lx

xkxk

.

L’equazione dei residui per un generico elemento di estremi 1x e nx è data da

( ) 0~~1

2

2

=

∂∂−

∂∂

∫ dx x xuu

i

x

xn

ψτ

, ni ,...,1= (3.9)

dove u~ (che verrà specificata tra breve) indica una generica soluzione

approssimata del problema. Dal momento che nella precedente equazione

compare un termine differenziale di secondo grado, appare opportuno intervenire

con un’integrazione per parti 4 del medesimo, trasformando la (3.9) in

nnnx

x

x

xi

ix

xi x

udx

dxd

xu

dxu

111

~~~∫∫

∂∂=

∂∂+

∂∂ ψ

ψψ

τ. (3.10)

Diversamente dai casi proposti nei paragrafi precedenti, compare, nel presente

problema, l’incognita temporale τ : la soluzione approssimata è cioè del tipo

( ) ( )∑=

=n

iii xaaxu

1

,;,~ τψτ . In questo modo, tuttavia, le iψ sono funzioni di due

3 Si veda il paragrafo 1.4.

4 dxdx

dxu

xu

dxxu nn n x

x

ix

x

x

xii ∫∫ ∂

∂−

∂∂=

∂∂

11 1

~~~2

2 ψψψ .

60

variabili, cui chiaramente non è possibile applicare gli elementi finiti

(monodimensionali) fin qui presentati. Per ovviare al problema, si opta per una

separazione di variabili che permette di considerare soluzioni del tipo

( ) ( ) ( )∑=

=n

iii xaaxu

1

;,~ ψττ .

A seguito di questa trasformazione le iψ sono tornate della forma consueta, ma

gli ia sono ora funzioni della variabile temporale τ . La soluzione sarà calcolata

in due passaggi distinti: il primo risolverà il problema rispetto alla variabile

spaziale, mentre il secondo rispetto a quella temporale.

Le derivate parziali di u~ , ∑=

=∂∂ n

jj

j

d

dau

1

~ψ

ττe ∑

=

=∂∂ n

j

jj dx

da

xu

1

~ ψ, sono sostituite in

(3.10), ottenendo

nnn

n n n

x

xij

n

j

x

x

ijjn

j

x

xij

x

x

n

j

x

x

x

x

n

ji

ijjij

j

xu

a dxdx

ddx

d

d

dadx

xu

dxdx

ddx

dadx

d

da

111

1 1 1

~

~

11

1 1

∂∂=

+

∂∂=+

∑ ∫∑ ∫

∫∑ ∫∑

==

= =

ψψψ

τψψ

ψψψψψ

τ

(3.11)

Il termine a destra dell’uguale non è stato trasformato, dal momento che è

esattamente l’espressione per le condizion i al contorno di Neumann: se esse

fossero presenti, sarebbe così possibile effettuare una sostituzione diretta. Nel

caso considerato, le condizioni al contorno sono formulate secondo Dirichlet, per

cui l’intero membro di destra può, in pratica, essere igno rato: verrà mantenuto

solo per correttezza formale. La (3.11) può essere riscritta, in modo più leggibile,

in forma di matrice:

61

( ) ( ) ( ) ( )

( ) ( ) ( ) ( )

( )

( )

( ) ( ) ( ) ( )

( ) ( ) ( ) ( )

( )

( )

∂∂

∂∂

=

+

+

∫∫

∫∫

∫∫

∫∫

n

n

nn

nn

nn

nn

x

xn

x

x

x

x

nnx

x

n

x

x

nx

x

nx

xnn

x

xn

x

xn

x

x

xu

xu

a

a

dxdx

xddx

xddx

dxxd

dxxd

dxdx

xddx

xddx

dxxd

dxxd

dda

dda

dxx xdxx x

dxx xdxx x

1

1

11

11

11

11

~

~1

2

1

1

111

1

1

111

ψ

ψ

τ

τ

ψψψψ

ψψψψ

ττ

ττ

ψψψψ

ψψψψ

(3.12)

Un analogo sistema di equazioni differenziali è ottenuto per ognuno degli N

elementi in cui il dominio è stato partizionato. La prima matrice è detta di massa,

mentre la seconda di rigidezza; essendo riferite al singolo elemento, sono dette

locali. Il vettore a destra dell’uguale è detto di flusso.

Scegliendo la funzione 0C -quadratica 2321

~ xxu ααα ++= , si ottengono, con

riferimento al primo elemento del dominio, le tre funzioni che lo descrivono 5:

5 Si veda il paragrafo 3.2.

62

( )( )( )( )

( )( )( )( )

( )( )( )( )2313

213

3212

312

3121

321

xxxxxxxx

xxxxxxxx

xxxxxxxx

−−−−

=

−−−−

=

−−−−

=

ψ

ψ

ψ

(3.13)

D’ora innanzi si farà riferimento e sclusivamente ad elementi 0C -quadratici e per

semplicità espositiva la matrice di massa dell’elemento k -esimo sarà indicata

come

( )

( ) ( ) ( )

( ) ( ) ( )

( ) ( ) ( )

=

kkk

kkk

kkk

k

mmmmmm

mmm

987

654

321

M ,

analogamente la matrice di rigidezza come

( )

( ) ( ) ( )

( ) ( ) ( )

( ) ( ) ( )

=

kkk

kkk

kkk

k

rrrrrr

rrr

987

654

321

R

e il vettore di flusso come

( )

( )

( )

( )

=

k

k

k

k

f

f

f

3

2

1

f .

63

3.3.2 Assemblaggio della matrice globale.

Il passo successivo consiste nell’assemblare le matrici locali in una matrice

globale, che rappresenti cioè il problema su tutto il dominio. Intuitivamente, per

ogni istante temporale, la soluzione globale è data dall’unione delle soluzioni

locali, essendo queste relative ad intervalli contigui del dominio. L’intuizione

tuttavia non basta, e, perchè il meccanismo di “giustapposizione” funzioni

correttamente, è necessario richiedere che, su ogni elemento, le funzioni iψ

rispettino due condizioni fondamentali:

• Condizione di completezza. La soluzione u~ e tutte le sue derivate, fino

all’ordine m (dove m2 è l’ordine dell’equazione differenziale del

problema originale), devono poter assumere qualsiasi valore finito

all’interno dell’elemento, al tendere a zero della dimensione di

quest’ultimo. Nel caso in cui la soluzion e abbia forma polinomiale di

grado p , ciò significa richiederne la completezza almeno fino al grado

m ( mp ≥ ).

• Condizione di continuità. In corrispondenza dei nodi di contatto tra gli

elementi, le funzioni u~ devono essere di classe 1−mC .

La prima condizione esercita la propria influenza nell’interno del singolo

elemento e garantisce la convergenza della soluzione approssimata alla soluzione

esatta. La seconda riguarda i nodi di contatto tra gli elementi e garantisce che, in

loro corrispondenza, non vi siano discontinuità tali da introdurre nuovi errori in

aggiunta a quelli “fisiologici” di discretizzazione. Le funzioni della famiglia (3.8)

rispettano entrambe le condizioni.

Nell’assemblaggio delle matrici locali, è necessario porre attenzione ai nodi di

contatto tra gli elementi: nel grafico seguente è riportata la corretta sequenza di

assemblaggio per i primi N elementi 0C -quadratici di un problema generico:

64

La sequenza kX , 11 +≤≤ Nk , indica i soli nodi esterni (riferiti agli N elementi),

mentre la sequenza ix , ( ) 1121 −+≤≤ Ni , considera tutti i nodi presenti. Per i

nodi di bordo vale la relazione 12 −= kk xX , mentre i nodi centrali rimangono

disaccoppiati.

Il relativo sistema assemblato in forma matriciale è rappresentato, per i primi due

elementi, da

2a

3a 4a

5a Na2

1a ( ) 112 −+Na

3x 2x 4x 5x 12 −Nx ( ) 112 −+Nx

1x Nx2

1X 2X 3X NX 1+NX

12 −Na

65

( ) ( ) ( )

( ) ( ) ( )

( ) ( ) ( ) ( ) ( ) ( )

( ) ( ) ( )

( ) ( ) ( )

( )

( )

( )

( )

( )

( ) ( ) ( )

( ) ( ) ( )

( ) ( ) ( ) ( ) ( ) ( )

( ) ( ) ( )

( ) ( ) ( )

( )

( )

( )

( )

( )

( )

( )

( ) ( )

( )

( )

fR

M

+=

++

+

+

23

22

21

13

12

11

5

4

3

2

1

29

28

27

26

25

24

23

22

21

19

18

17

16

15

13

12

5

4

3

2

1

29

28

27

26

25

24

23

22

21

19

18

17

16

15

14

13

12

11

00

00

00

00

00

00

00

00

f

f

ff

f

f

a

a

a

a

a

rrr

rrr

rrrrrr

rrr

rrr

dda

dda

dda

dda

dda

mmm

mmm

mmmmmm

mmm

mmm

14

11

τ

τ

τ

τ

τ

ττ

ττ

ττ

ττ

ττ

(3.14)

o, in forma compatta,

( ) ( ) fRaaM =+ τττ

dd (3.15)

La soluzione assemblata, per ciascun istante di tempo, è allora data

( ) ( ) ( ) ( )( ) ( ) ( ) ( )( )

( ) ( ) ( )( ) ( ) ( ) ( )( )

elemento ultimo

NN

NN

NN

esimo-e elemento

ei

ei

ei

elemento secondoelemento primo

aaaaaa...

aaaaaaaxU

3121131211

235

224

213

133

122

111

...

...;,~

ψψψ

ψψψψψψτ

+−+− +++++++

++++++=

(3.16)

66

dove ( )eiψ , 3,2,1=i , indica la funzione iψ , rappresentata nella relazione (3.13),

calcolata in riferimento all’ e -esimo elemento.

Se il problema originale fosse indipendente dal tempo, il sistema (3.14) sarebbe

algebrico e non differenziale: la sua soluzione fornirebbe direttamente i valori dei

parametri ia da sostituire in (3.16) per calcolare la soluzione.

Per ottenere la soluzione finale del problema è necessario dunque risolvere il

sistema differenziale rispetto al tempo (3.14): per farlo ci si può avvalere di uno

schema alle differenze finite. Optando per uno schema backward del tipo

presentato nel paragrafo 2.2.2, la (3.15), riferita al generico istante τ , diviene:

( ) ( )( ) ( ) fRaaaM =++∆

−+1

1 ττ

ττ (3.17)

la cui soluzione è

( ) ( )( )fMaRMa +

∆+=+

−

ττ

τ1

11 . (3.18)

Il vettore ( )τa è deputato a introdurre le condizioni dell’equazione differenziale

di partenza: i termini del vettore iniziale ( )0a sono costituiti da

( ) ( ) ( )( )0,max0 1121

21

ii xkxki ee a +− −=

mentre i termini estremi del vettore rappresentano le condizioni al bordo

( ) 0,lim =±→

τxuLx

:

( ) ( ) 01 == ττ naa .

La sostituzione in (3.16) dei valori ia permette di trovare le soluzioni cercate.