Indice Prefazione 1 CLASSIFICAZIONE DELLE EQUAZIONI ALLE DERIVATE PARZIALI.............................................................................4 1.1 Equazioni alle derivate parziali del secondo ordine.................................4 1.2 Sistemi del primo ordine in due variabili................................................5 2 METODI VARIAZIONALI........................................................................9 2.1 Elementi di calcolo delle variazioni........................................................9 2.2 Formulazione debole e formulazione variazionale di problemi ellittici....................................................................................11 2.3 Principi variazionali nei problemi agli autovalori.............................................................................................15 2.4 Metodi variazionali approssimati........................................................16 2.4.1 Metodo di Rayleigh-Ritz.............................................................16 2.4.2 Metodo di Galerkin.....................................................................17 2.4.3 Metodo di Rayleigh-Ritz nei problemi agli autovalori.....................................................................................18 3 ESEMPI APPLICATIVI..........................................................................20 3.1 Equazione della membrana.................................................................20 3.1.1 Deformazione della membrana sotto l’azione di un carico ripartito....................................................................20 3.1.2 Autovibrazioni della membrana..................................................23 3.2 Svuotamento a potenziale di un serbatoio...........................................27 3.3 Oscillazioni trasversali di una trave.....................................................33 3.3.1 Metodo di Rayleigh-Ritz per la trave.........................................38 3.4 Equazione della piastra.......................................................................39 3.4.1 Deformazione della piastra sotto l’azione di un carico ripartito.....................................................................39 3.4.2 Autovibrazioni della piastra.........................................................41 4 CENNI SUL METODO AGLI ELEMENTI FINITI (FEM).....................45 4.1 Considerazioni generali........................................................................45 4.2 Costruzione generale della procedura numerica...................................45 4.3 Caso 2D con mesh ad elementi triangolari e piecewise lineari.............46 5 METODI ALLE DIFFERENZE FINITE.................................................50 5.1 Considerazioni generali......................................................................50 5.2 Analisi degli schemi numerici..............................................................51 5.2.1 Consistenza.................................................................................51 5.2.2 Convergenza................................................................................53 5.2.3 Stabilità......................................................................................54 5.3 Schemi espliciti e schemi impliciti......................................................56 5.4 Condizione di stabilità di Von Neumann............................................56 5.5 Dissipazione e Dispersione.................................................................58 6 APPENDICE: RICHIAMI DI MATEMATICA E APPROFONDIMENTI........................60 6.1 Formule di Green-Gauss......................................................................60 1
Il presente testo raccoglie le lezioni da me tenute per il corso di Meccanica
Computazionale nell’ambito della Laurea Magistrale in Ingegneria Civile per la
Protezione dai Rischi Naturali. Il settore scientifico disciplinare in cui si inserisce
tale corso è MAT/07 (Fisica Matematica).
L’obiettivo di tale insegnamento è quello di fornire allo studente della Laurea
Magistrale gli elementi di base per inquadrare e risolvere, dal punto di vista
teorico e computazionale, alcuni problemi fisico-matematici che si incontrano
sia nell’ingegneria strutturale che in quella idraulica.
Le tecniche di risoluzione numerica utilizzate sono: i metodi variazionali
approssimati, il metodo agli elementi finiti (FEM) e i metodi alle differenze
finite.
La scelta del contesto applicativo si focalizza su problemi che, sebbene sem-
plici dal punto di vista della geometria che li descrive, permettano di eviden-
ziare aspetti dinamici fondamentali per un corretto inquadramento concettuale
di problemi più complessi. Essendo il corso destinato a studenti di ingegneria,
molti aspetti matematici sono stati volutamente semplificati, a volte a scapito
del rigore.
Infine tutti i problemi applicativi sono stati risolti utilizzando il programma
di calcolo simbolico e numerico Mathematica. Tale strumento permette una
implementazione rapida e compatta dei programmi e, tramite le efficacissime
animazioni, fornisce uno strumento didattico versatile e stimolante per la for-
mazione di una corretta immagine mentale dei problemi affrontati.
Roma 2015
Giampiero Sciortino
3
1 Classificazione delle equazioni differenziali alle
derivate parziali
1.1 Equazioni alle derivate parziali del secondo ordine
Considereremo equazioni lineari o quasi lineari, in cui i coefficienti dell’equazione
possono dipendere dalla funzione incognita (ma non dalle sue derivate) e dalle
variabili indipendenti. La struttura della generica equazione nella incognita:
= (x) ∈ x+(1 2 ) ∈ sarà:
X=1
2
+
X=1
+ = (1)
Stante la 2
= 2
si può supporre senza scapito di generalità che la
matrice = (()) dei coefficienti dei termini del secondo ordine sia simmetrica.
I suoi autovalori saranno allora tutti reali. Denotiamo con + − 0 rispetti-vamente il numero di quelli positivi, negativi e nulli in un certo x = x0 = 0
Ovviamente = ++−+0 Diremo allora che in x = x0 = 0 l’equazione
(1) è di tipo (per brevità con ∧∨ indicheremo rispettivamente la congiunzione”e” e la disgiunzione ”oppure”):
Si chiama superficie caratteristica una superficie di equazione (x) =0 che
verifica la relazione
(∇z∇z) =X
=1
= 0 (2)
dove (• •) indica il prodotto scalare.Nel caso in cui = 2, ossia il problema sia in due variabili la equazione (1)
è di tipo:
2
2+ 2
2
+
2
2+ termini10 = (3)
Se alla equazione (3) associamo il problema di Cauchy, consistente nell’assegnare
i valori di e della sua derivata in direzione normale (n) ad un’assegnata
curva Γ, si può dimostrare che tale problema è indeterminato solo su delle linee
eccezionali dette linee caratteristiche che soddisfano l’equazione differenziale
ordinaria ottenibile dalla (2)
( )
µ
¶µ
¶= 0
2 + 2 + 2 = 0
(4)
4
Passando da una relazione implicita ( ) = 0 ad una esplicita = ()⇒( ()) = 0 ⇒ − = ed inserendo nelle (4), segue l’equazione
differenziale ordinaria delle linee caratteristiche:
()2 − 2() + = 0 (5)
Le precedenti forniranno una famiglia di due curve reali e distinte, coinci-
denti, complesse coniugate, rispettivamente se 2 − T 0 Poiché la matrice associata alla equazione (3) è:
=
µ
¶segue che:
det(− ) = 2 − +∆ = 0
12 =±√2−4∆
2
+ + ∆ + − 2 2 − 4∆ = (− )2 + 2 ≥ 0
Supposto ad es. 0 segue che se ∆ 0 ⇔ 2 − 0 ⇒ 12 0
ed il sistema è Ellittico (entrambi gli autovalori sono positivi). Se ∆ = 0
segue che un è nullo e il sistema è Parabolico. Infine se ∆ 0 un auto-
valore è positivo e uno negativo e il sistema è Iperbolico. Se = 0 ⇒ =
− ⇒ ∆ = 2 − = 2 + 2 0 ed il sistema è ellittico. (Analoghe consider-
azioni se 0) Da quanto mostrato, segue che le equazioni Ellittiche in due
variabili sono caratterizzate da una famiglia di linee caratteristiche complesse
coniugate, le Iperboliche da una famiglia di linee reali e distinte, le Paraboliche
da una famiglia di linee reali e coincidenti. Le equazioni ellittiche modellano
in genere processi stazionari. Le Iperboliche traducono in genere principi di
conservazione con reversibilità temporale. Le paraboliche processi che evolvono
nel tempo in modo unidirezionale, nel senso della irreversibilità. Una inversione
della direzione del tempo porta a problemi mal posti.
1.2 Sistemi del primo ordine in due variabili
Consideriamo il sistema di equazioni in due variabili indipendenti :
U
+
U
= C (6)
dove:
U = [1( ) 2( ) ( )]
= (()) = (())
C = [1 2 ]
= 1 2
5
La classificazione del sistema (6) si effettua in relazione alle proprietà del seguente
problema generalizzato agli autovalori ( + [12 ] ):
=
() + det(− ) (7)
supponendo che almeno una delle matrici sia non singolare. Suppor-
remo nel seguito che det 6= 0.Diremo allora che il sistema (6) è:
ellittico se gli zeri di () = 0 sono tutti complessi
iperbolico se gli zeri di () = 0 sono tutti reali e il problema agli auto-
valori (7) ammette una base di autovettori
parabolico se gli zeri di () = 0 sono tutti reali e il problema agli auto-
valori (7) non ammette una base di autovettori.
Se le matrici e e il vettore C dipendono da e il sistema (6) è
lineare. Se dipendono anche da U il sistema è detto quasi lineare. Nel primo
caso la classificazione può essere punto-dipendente nel secondo anche soluzione-
dipendente.
Nel caso di sistema iperbolico, (denoteremo in seguito con la variabile
(variabile tempo)) la (7) diventa per ipotesi di non singolarità di :
−1 = (8)
cosicchè nella base degli autovettori la matrice = −1 è diagonale:
e = Λ (9)
dove Λ = (1 2 ) è la matrice diagonale che contiene gli auto-
valori.
Posto allora:
= −1 = −1 e (10)
dove è una matrice non-singolare che fa passare dalla base degli autovettori
alla base canonica cui si riferiscono le componenti di. Il significato invece della
matrice −1 è immediato tenendo conto che, la sua −esima colonna, contienele componenti dell’−esimo autovettore del problema (8).Il sistema (6) moltiplicando a sinistra per −1 diventa:
−1U
+
U
= −1C (11)
e quindi per le (10):
−1 eU
+
U
= −1C⇒ ΛU
+
U
= −1C (12)
Consideriamo adesso alcuni casi particolari.
6
Se le matrici e sono costanti allora lo è anche e possiamo effettuare
un cambio della variabile dipendente U ponendo:
eU = U⇒ eU( )
= U
( )
cosicchè le (12) diventano:
eU
+ Λ eU
= −1C (13)
o anche scritte per componenti:
f
+ e
= () e = (−1C)
dove ()(·) + (·) + (·) rappresenta la derivata lungo la − linea caratteristica di pendenza = , eU =
he1 e2 ei edessendo (−1C) la − componente del vettore −1C.Qualora nelle precedenti C = 0 segue che le f non variano lungo le cor-
rispondenti linee caratteristiche di pendenza = , e il sistema può essere
integrato analiticamente ottenendo:
e = (− ) (14)
= 1 2
con funzioni arbitrarie di classe 1
Le precedenti ci dicono che le e( ) non variano lungo le linee − =cost,che sono le linee caratteristiche del sistema e le e( ) costituiscono quindi degliinvarianti.
Talvolta le matrici e non sono costanti ma dipendono unicamente da
U, come in alcuni problemi quasi lineari. In tal caso, se sussiste la possibilità
di porre (tale possibilità non è vera in generale):
U
=
F(U)
U
=
F(U)
⇔ =
F(U)
U(15)
F(U) + [1[U]2[U] [U]]
le (12) diventano:
F(U)
+ Λ
F(U)
= −1C (16)
o anche scritte per componenti:
(U)
+
(U)
= () = (−1C)
7
Anche in tal caso se C = 0 otteniamo che le [U] non variano lungo le
linee caratteristiche di equazione = e quindi sono degli invarianti detti
invarianti di Riemann. In tal caso però gli autovalori dipendono daU e quindi
le linee caratteristiche non sono tracciabili a priori. La possibilità comunque di
costruire invarianti è legata alla possibilità di soddisfare le (15), ossia che la
matrice sia lo Jacobiano di qualche F(U).
Se comunque una certa equazione del sistema (12) con C = 0, diciamo l’− equazione:
X=1
(+ ) =
X=1
() = 0
può essere scritta come:
[U]
+
[U]
= () = 0
allora [U] è un invariante associato alla caratteristica di pendenza =
Nel caso più generale, il sistema (12) può essere semplicemente scritto facendo
intervenire le derivate delle lungo le linee caratteristiche nel seguente modo:
X=1
() = (−1C)
= 1 2
8
2 Metodi Variazionali
2.1 Elementi di calcolo delle variazioni
Sia Ω ⊂ un insieme aperto (campo). In Ω definiamo il seguente spazio di
funzioni, detto spazio delle funzioni a quadrato sommabile in Ω così definito:
2(Ω) =
½()
¯ ∈ () ∈
ZΩ
()2Ω ∞¾
Tale spazio è uno spazio vettoriale sul campo dei reali, nel quale possiamo
definire il seguente prodotto scalare fra due suoi elementi:
(() ()) +ZΩ
()()Ω
Con Ω = Ω ∪ Ω indicheremo la chiusura di Ω ossia l’unione con la sua
frontiera.
Una applicazione , che ad ogni funzione appartenente ad un assegnato
spazio A, associa un numero reale [], si chiama funzionale. In simboli:
: A→
()→ []
In molti problemi, A ⊂2(Ω), ossia A è un sottoinsieme dello spazio vetto-
riale 2(Ω) Tale applicazione è un funzionale di una sola variabile Qualora
l’applicazione sia del tipo:
: A1 ×A2 × ×A →
(1() 2() ())→ [1 2 ]
dove () ∈ A = 1 2 si parla di funzionale di più variabili.
Per definire correttamente molti problemi applicativi è utile introdurre i
seguenti sottospazi vettoriali di 2(Ω) :
(Ω) = l’insieme delle funzioni di 2(Ω) che hanno derivate continue in Ω
fino all’ordine = 0 1 2
(Ω) =l’insieme delle funzioni di 2(Ω) che hanno derivate fino all’ordine
= 1 2 3appartenenti ancora a 2(Ω)
Andiamo ora a definire il concetto di variazione di un funzionale. In sostanza
vogliamo calcolare la variazione che subisce il funzionale quando si passa da una
funzione ad un’altra funzione +∆ ossia [+∆]− [] Per motivi che inseguito saranno chiari, conviene scrivere ∆ = dove è uno scalare e una
funzione. Supponendo come abbiamo detto che il funzionale sia definito in Anon è affatto detto che se ∈ A allora + ∈ A (questo avviene solo se Aè uno spazio vettoriale). Per essere allora sicuri che + ∈ A, introduciamoun insiemeM (che si può dimostrare essere uno spazio vettoriale) così definito:
M = () |+ ∈ A∀ ∈ A∀ ∈
9
Come vedremo in alcuni problemi, quando A è un sottospazio vettoriale di
2(Ω) allora si può sempre supporre M = A, mentre in altri problemi i dueinsiemi sono distinti.
Se esiste il seguente limite:
lim→0
[+ ]− []
=
¯=0
[+ ]∀ ∈M (17)
chiameremo tale limite la variazione prima di nel punto in direzione e
la indicheremo con la notazione [; ] Analogamente a quanto avviene per le
funzioni, si può dimostrare che, se in 0 il funzionale ha un massimo o minimo
relativo e se in 0 ammette variazione prima, allora:
[0; ] = 0∀ ∈M (18)
Nel caso di funzionale di più variabili, si definisce una variazione prima
parziale rispetto ad ogni variabile = 1 2 ponendo:
[1 2 ; ] +
lim→0
[1 2 + ]− [1 2 ]
∀ ∈M
Nel caso di funzioni ordinarie, sussiste la proprietà che la derivata lungo una
assegnata direzione è data dal prodotto scalare del suo gradiente per il versore in
quella direzione. Anche per i funzionali sussiste una proprietà analoga. Si può
infatti dimostrare che, sotto determinate condizioni, [; ] che altro non è che
la derivata di nella direzione , può essere scritto come [; ] = ( ) dove
è una funzione (dipendente da ) appartenente a 2(Ω) Per quanto detto è
naturale definire tale funzione il gradiente di e scrivere
[; ] = (∇ [] ) (19)
Ora, per la (18), se in 0 il funzionale ha un massimo o minimo relativo
si ha: (∇ [0] ) = 0∀ ∈M ⇒ ∇ [0] = 0 Tale equazione prende il nomedi equazione di Eulero associata al funzionale [] Quando, come vedremo nel
prossimo capitolo, ad una equazione alle derivate parziali si associa un equiv-
alente problema variazionale, in sostanza si definisce il funzionale in modo che
l’equazione di Eulero ad esso associata sia l’equazione alle derivate parziali di
partenza.
Un caso particolarmente semplice di funzionale con diverse applicazioni in-
teressanti, è quello definito dalla relazione:
[] =
Z 2
1
( () 0()) (20)
dove ( 0) è una assegnata funzione e 1,2 sono assegnati. Il problemaconsiste nel trovare una funzione () ∈ A =1 [1 2] che minimizzi o mas-simizzi [] e che passi per i punti (1 1),(2 2). Imponendo l’annullamento
ossia il più piccolo autovalore 1 è il valore minimo che assume il funzionale
[] al variare di in A0 escludendo la funzione identicamente nulla e il valoreche assume [] sulla generica autofunzione del problema (40) fornisce il
corrispondente autovalore come si verifica immediatamente moltiplicando
per l’equazione del problema (40) e integrando su Ω. Si verifica inoltre che
[] =2[]
2= [] (ossia il funzionale è omogeneo di grado zero). Tale
proprietà garantisce da un lato che che la ricerca del minimo di [] in A0−0è equivalente alla ricerca del minimo del numeratore della (42) col vincolo che il
denominatore valga uno (vedi Appendice), dall’altro che la relazione = []
non presenti ambiguità, dato che le autofunzioni sono definite a meno di un
fattore moltiplicativo.
Infine ogni ( = 2 3 ) costituisce un minimo vincolato per [] nel
senso che:
= min∈(A0−0)∩M
[] (44)
M +n ∈ 2(Ω)
¯(12 12) = 0 = 1 2 − 1
o(45)
ossia è il minimo di [] in A0 − 0 col vincolo che sia ortogonale(rispetto al peso ) alle autofunzioni =12−1.
2.4 Metodi variazionali approssimati
2.4.1 Metodo di Rayleigh-Ritz
Tale metodo si applica a quei problemi che ammettono una formulazione vari-
azionale, ossia problemi in cui si cerca una funzione che minimizzi (o mas-
simizzi o renda stazionario) un certo funzionale [] in un definito spazio AD’ora in poi, senza specificare massimo, minimo o punto stazionario, diremo
semplicemente punto estremale intendendo un punto che annulla la variazione
prima del funzionale. Consideriamo allora il caso del problema (24) che è ab-
bastanza generale, e gli spazi definiti dalle (28) che per comodità riscriviamo (al
ªSia ora noto il funzionale []. Sia poi 0 ∈ A una arbitraria funzione
assegnata (se è nota anche in Ω o se è prolungabile in Ω conviene assumere
0 = ) e siano 1 2 funzioni linearmente indipendenti appartenenti
aM. Allora la funzione:
() = 0() +
X=1
() (47)
appartiene ad A ∀(1 2 ) ∈
Al variare comunque di (1 2 ) ∈ la funzione spazzerà un
insieme A ⊂ A. Il metodo di Rayleigh-Ritz consiste nell’approssimare la
ricerca dell’estremale di in A, nella ricerca dell’estremale di in A Tanto
meglio lo spazio A spazzato dalla approssima bene A, tanto più precisoè il risultato fornito da tale metodo. Ovviamente trovare un estremale di in
A equivale a trovare i punti estremali della seguente funzione delle variabili
(1 2 ) ∈ :
Λ(1 2 ) + [0 +
X=1
] = []∈A
(48)
ossia risolvere il sistema di equazioni in incognite:
Λ = 0 = 1 2 (49)
Trovati i coefficienti che soddisfano le (49), trovo () che è la funzione
approssimante la soluzione esatta del problema variazionale.
In molti casi di interesse applicativo, data la struttura quadratica dei fun-
zionali, il sistema (49) è lineare e quindi facilmente risolvibile dal punto di vista
computazionale.
2.4.2 Metodo di Galerkin
Tale procedura si utilizza per approssimare soluzioni deboli, quando la formu-
lazione variazionale non è possibile. Il problema è allora trovare una funzione
∗ ∈ A tale che
[∗ ] = []∀ ∈M (50)
Siano allora 1 2 , funzioni linearmente indipendenti, appartenenti
aM, per la costruzione di una approssimante ∗ sempre secondo una formula
analoga alla (47):
∗ () = 0() +
X=1
∗() (51)
17
e siano 1 2 funzioni linearmente indipendenti appartenenti a
M (normalmente si assumono coincidenti con le 1 2 ) necessarie ad
approssimare la condizione che impone il soddisfacimento della [∗ ] = []
per ogni ∈M. Scrivendo allora:
[∗ ] = [] = 1 2 (52)
otteniamo il sistema di equazioni in incognite ∗1 ∗2
∗ Imporre
che la (50) sia soddisfatta per = 1 2 è sufficiente affinchè valga per
tutte le ∈M =spazio vettoriale a dimensioni spazzato dalle 1 2 .
Risolto il sistema, dalla (51) ricaviamo ∗ () che è la funzione approssimantela soluzione esatta del problema debole (50). Si può poi dimostrare che tale
procedura coincide con quella di Rayleigh-Ritz, qualora la formulazione debole
(50) ammette una formulazione variazionale equivalente.
2.4.3 Metodo di Rayleigh-Ritz nei problemi agli autovalori
Applichiamo la procedura di Rayleigh-Ritz al problema agli autovalori (40).
Siano 1 2 , funzioni linearmente indipendenti, appartenenti allo spazio
azionali approssimati, la funzione approssimante si deve porre nella forma:
2 = 0 +
X=1
(86)
0 ∈ A ∈M
La scelta delle , tenendo conto che l’appartenenza ad M richiede che tali
funzioni si annullino sulla frontiera Ω può essere fatta utilizzando le funzioni:
= sin(
) sin(
) (87)
Per quanto riguarda la scelta di 0tale scelta è arbitraria ma l’appartenenza
ad A richiede che 0 assuma su Ω i valori di E’ allora naturale cercare di
vedere se è prolungabile in Ω Ora sappiamo che 1() è il valore che assume
per = , mentre 3() il valore che assume per = 0 Per prolungarla in
maniera continua all’interno del dominio fluido, è naturale provare a connettere
tali due funzioni attraverso una combinazione lineare a coefficienti dipendenti
da ponendo:
= 0 =(− )
3() +
1() (88)
Come è immediato verificare tale funzione è definita e continua in tutto
Ω ∪ Ω e su Ω assume i valori di Va sottolineato che tale funzione è solo
continua, ma ha derivate discontinue come . L’appartenenza ad A non richiedeperò alcuna condizione di continuità delle derivate ma solo che queste siano a
quadrato integrabile (appartenenza ad 1(Ω)), condizione senz’altro verificata
data la finitezza delle derivate nei punti dove esistono.
Inseriamo ora l’approssimante (86) nel funzionale (85) e cerchiamo gli estremi
della funzione:
Λ = [0 +
X=1
] (89)
al variare di ∈ Imponendo l’annullamento delle derivate parziali
Λ, segue, dopo alcuni passaggi, che i coefficienti devono verificare
il sistema lineare:
X=1
¿
+
À= −
¿0
+
0
À(90)
= 1 2
29
L’ortogonalità delle in Ω ha come conseguenza che nella sommatoria a
primo membro delle (90) solo i termini con = e = sono non nulli cosicchè
il sistema è immediatamente risolvibile in quanto diagonalizzato e porge:
=
+ −¿0
+
0
À(91)
+¿
+
ÀEseguendo i calcoli si trova che:
= −
(− )(sin(
)− sin(
))
=2
4(2
+
2
)
Nelle Fig.9 e Fig.10 sono riportati i grafici delle linee di livello (linee di
corrente) per l’approssimante (86) con = 20, relativi rispettivamente ad un
serbatoio quadrato di lato unitario con = 1 e con efflusso compreso nella
zona fra 04 06 e 06 08:
Fig. 9 Fig.10
Prima di concludere tale esempio facciamo una interessante osservazione.
L’approssimante (86), avendo come addendo 0 è una funzione continua ma
con derivate discontinue. Per → ∞ deve però fornire la soluzione esatta
dell’equazione di Laplace relativa al problema (83), che richiede funzioni con
derivate continue fino al secondo ordine (appartenenza a 2(Ω)). Per visualiz-
zare come la funzione 2 "transmigra" per → ∞ da 1(Ω) verso 2(Ω)
30
(equivalenza fra la formulazione debole e quella classica) osserviamo i grafici di
Fig.11 e Fig.12 che riportano l’andamento di 2( 2) plottato al variare
)(1 H )(2 C
2N
di per = 2 rispettivamente per = 0 = 5 = 20 = 40 :
= 0 = 5
Fig11
= 20 = 40
Fig.12
Come si vede, all’aumentare del numero di modi, la funzione perde le dis-
continuità e tende a diventare liscia come richiede la formulazione classica del
problema (83).
31
Lo svuotamento a potenziale del serbatoio, può essere fatto anche in termini
della funzione potenziale risolvendo sempre l’equazione di Laplace ma con
condizioni al contorno alla Neumann essendo in tal caso V =∇. Al posto delproblema (83) avremo allora il problema:½ ∇2 = 0 in Ω
= in Ω(92)
dove deve valere (la normale n è assunta esterna):
=
⎧⎪⎪⎪⎪⎪⎪⎨⎪⎪⎪⎪⎪⎪⎩
= −1 in 1 = 0 ≤ ≤ = = 0 in 2 = 0 ≤ ≤ = − = 0 in 4 = 0 ≤ ≤ = 0
− =⎧⎨⎩ 0 per 0 ≤ ≤
3 per
0 per ≤ ≤
in 3
Da notare che tale funzione verifica la condizione di risolubilità (39) in
quanto: ZΩ
Ω = 3(− )− 1 = − = 0
Dalla formula (36) sappiamo che il problema (92) è equivalente a trovare
Poichè in tale problema non è richiesto nessun comportamento speciale
sulla frontiera Ω per appartenere ad A, l’approssimante conviene definirlaattraverso una serie di coseni che non si annullano identicamente sulla fron-
tiera. Se si scegliesse una serie di seni che si annullano su Ω come nel prob-
lema (83), il termine
ZΩ
nella (93) sarebbe identicamente nullo calcolato
sull’approssimante e il sistema lineare per calcolare i coefficienti dello sviluppo
sarebbe quindi omogeneo ammettendo la sola soluzione nulla. Posto allora:
2 =
X=0
Λ + [2 ] (94)
+ cos(
) cos(
) 00 + 0
imponendo l’annullamento delle derivate parziali Λ si ottiene un sis-
tema lineare nelle incognite , i cui valori forniscono l’approssimante (94).
L’aver posto 00 + 0 è conseguenza del fatto che l’equazione di Laplace, concondizioni alla Neumann, ammette infinite soluzioni che differiscono per una
32
costante arbitraria e quindi l’annullamento di 00 equivale ad aver fatto una
particolare scelta per tale costante arbitraria.
Nelle Fig.13 sono riportate le soluzioni sovrapposte alle analoghe figure 9 e
10 per visualizzare la rete del moto.
Fig.13
3.3 Oscillazioni trasversali di una trave
Consideriamo una trave omogenea di lunghezza , sezione , modulo elastico
e densità , vincolata ai due estremi mediante una fissata combinazione di due
dei tre seguenti vincoli: appoggio. incastro o estremo libero. Per descrivere le
oscillazioni forzate trasversali di piccola ampiezza della trave, bisogna risolvere
l’equazione del quarto ordine:
22 + 244 = ( ) (95)
( ) + ( )()
dove 2 + (), ( ) è il carico distribuito (positivo verso l’alto),
il momento d’inerzia rispetto all’asse neutro e ( ) la deformata istantanea
della trave rispetto alla posizione di riposo.
z
F
x
33
Cominciamo con studiare le autovibrazioni della trave, ossia oscillazioni non
forzate in cui ogni elemento della trave oscilla con una medesima frequenza.
Dall’equazione (95) annullando allora il carico distribuito e cercando soluzioni
del tipo:
( ) = () (96)
segue il problema agli autovalori del quarto ordine:⎧⎨⎩ 44 = , + 22
11 = 22 = 0 = 0
11 = 22 = 0 =
(97)
consistente nel trovare soluzioni non identicamente nulle assieme al cor-
rispondente autovalore 0.
Gli interi 1 2 1 2 assumono i valori 0 1 2 3 (00 è da intendersi
come ) in relazione alle condizioni di vincolo in = 0 e = . In partico-
lare per l’appoggio si devono annullare lo spostamento (indice 0) e il momento
(indice 2), per l’incastro spostamento e rotazione (indici 0 e 1) e per l’estremo
libero momento e taglio (indici 2 e 3). Il problema (97) può essere risolto analiti-
camente nel seguente modo. La soluzione generale della equazione differenziale
ordinaria del quarto ordine 44 = si verifica immediatamente essere:
= sin() + sinh() + cos() + cosh() (98)
dove + 4√ e sono quattro costanti arbitrarie. Imponendo
le quattro condizioni al contorno che impongono il soddisfacimento dei vin-
coli, si ottiene un sistema omogeneo di quattro equazioni lineari nelle incog-
nite . Annullando il determinante di tale sistema, al fine di im-
porre l’esistenza di autosoluzioni (non identicamente nulle) si perviene ad una
equazione del tipo:
() = 0 (99)
le cui infinite soluzioni reali e positive 1 2 forniscono le radiciquarte degli infiniti autovalori del problema (97). Per ogni autovalore è possibile
ricavare risolvendo il sistema omogeneo le (a meno di una costante
moltiplicativa) che inserite nella (98), forniscono la corrispondente autofunzione,
detta anche modo di oscillazione. Autovalori e autofunzioni possono essere
ordinate per valori crescenti degli autovalori:
(1 (1)
) (2 (2)
) ( ()
) (100)
1 2 3
poiché si può dimostrare che esiste un autovalore più piccolo 1.
Ad ogni autovalore corrisponde la pulsazione = √ = 2, che è la
−esima pulsazione di risonanza della trave. La prima frequenza di risonanza èdunque la più bassa, mentre le altre sono via via crescenti con l’ordine della riso-
nanza stessa. Le funzioni ()
() sono particolari soluzioni (particolari in
34
quanto ogni elemento della trave oscilla con uguale pulsazione ) dell’equazione
omogenea che si ottiene dalla (95) ponendo = 0, ossia dell’equazione:
22 + 244 = 0 (101)
L’insieme di tutte le soluzioni di tale equazione, essendo questa lineare ed
omogenea, costituisce uno spazio vettoriale S come facilmente si verifica. Fra
due elementi (funzioni) 1 2 ∈ S (ossia soddisfacenti la (101)), definiamo il
seguente prodotto scalare:
(1 2) +Z0
12 (102)
Se (1 2) = 0 diremo che le due funzioni sono ortogonali, in analogia ai
vettori dello spazio ordinario.
Dimostriamo ora la seguente proprietà fondamentale. Per le ipotizzate con-
dizioni di vincolo alle due estremità della trave (appoggio, incastro o estremo
libero) i modi di oscillazione sono fra loro ortogonali, ossia:
(()
()
) = 0 se 6= (103)
Ovviamente poiché si sono scartate le soluzioni banali (ossia identicamente
nulle), sarà sempre (()
()
) 0∀ = 1 2 3.Per dimostrare la (103),
cominciamo con l’osservare che ∀1 2 ∈ S ⇒ (414 2) = (1
424).
Tale relazione si ricava facilmente integrando quattro volte per parti fino a
“spostare” la derivata quarta sulla funzione 2. Ad ogni integrazione per parti,
tenendo conto delle possibili condizioni di vincolo ipotizzate, spariscono sempre
i termini variati fra 0 ed . Ricordando allora che una generica ()
verifica la
4()
4 = ()
, segue allora che:
(4()
4 ()
) = (()
()
) = (()
4()
4) = (()
()
)
e dunque (()
()
) = (()
()
) ⇒ ( − )(()
()
) = 0. Se
6= allora ( − ) 6= 0⇒ (()
()
) = 0 il che dimostra la (103).
L’insieme delle autofunzioni ortogonalin()
()o=123
può essere usata
ora quale base per rappresentare un generica funzione ( ) che deve soddisfare
l’equazione di partenza non omogenea (95). A tale scopo basta porre:
( ) =
∞X=1
()()
() (104)
essendo () dei coefficienti incogniti del tempo, che rappresentano l’ampiezza
del -esimo modo di oscillazione ()
(). Inserendo la (104) nella (95) segue:
∞X=1
(2()
2()
() + 2()4
()
()
4) = ( )
35
Ricordando al solito che ()
verifica la 4()
4 = ()
, con = 22,
segue allora la:∞X=1
(2()
2+ 2())
()
() = ( ) (105)
Moltiplicando scalarmente entrambi i membri della precedente per ()
e
ricordando la proprietà di ortogonalità (103), segue:
2()
2+ 2() = () (106)
() +f()
(()
()
) f() + (( ) () ())
Qualora il carico ( ) = ( ) sia concentrato in un punto (even-
tualmente mobile ()), le precedenti equazioni possono essere ancora utilizzate
attraverso un opportuno passaggio al limite. Sia infatti () un carico con-
centrato in (). Immaginiamolo ora distribuito uniformemente nell’intervallo
[()− 2 ()+ 2] con un carico distribuito uniforme di intensità ( ) =
( ) = (). Per → 0 tale carico distribuito tenderà al carico concen-
trato () in (), col che il termine (( ) ()
()) nel passaggio al limite
diventerà:
lim→0(( )
()
()) = lim→0
1
(()
()
()) =
=()
lim→0
1
()+2Z()−2
()
() =()
()
(())
avendo applicato il teorema della media nell’ultimo passaggio.
Segue che per un carico () concentrato in () le equazioni (106) diven-
tano:
2()
2+ 2() =
()
()
(())
(()
()
)
Generalizzando e sovrapponendo gli effetti per il caso di un carico distribuito
( ) e carichi concentrati () () in
() () ( = 1 2 ) avremo le
equazioni:
2()
2+ 2() = () (107)
() +f()
(()
()
) f() +
1
( ( )
()
()) +1
X=1
() ()()
(() ())
L’ampiezza () del -esimo modo di oscillazione ()
(), soddisfa dunque
l’equazione di un oscillatore armonico di pulsazione propria , eccitato dalla
36
forzante () proporzionale alla proiezione del carico sul modo stesso ()
().
Se tale proiezione è nulla, il modo corrispondente non viene eccitato. Se la ()
oscilla con pulsazione prossima alla , allora il -esimo modo di oscillazione
()
() va in risonanza. Se infine il carico non dipende esplicitamente dal tempo,
allora le () = sono delle costanti. Le (107) ammettono allora le soluzioni
costanti = 2 che inserite nella (104) forniscono la deformata della trave
(soluzione stazionaria) sotto l’azione statica del carico. Nel caso più generale,
l’integrazione delle (107) fornisce delle funzioni del tempo () che inserite
nella (104) forniscono la deformata istantanea della trave ( ). Risolvendo
ad esempio le (107) con condizioni iniziali (0) = 0 0(0) = 0, si sta impo-
nendo che all’istante zero la trave sia indeformata e ferma. Relativamente a tale
caso, supponiamo che il carico solleciti la trave, a partire dall’istante = 0, per
un intervallo di tempo finito. Valutiamo allora l’energia che tale sollecitazione
trasferisce alla trave e come tale energia si distribuisce tra i vari modi di oscil-
lazione. L’energia della trave è la somma di quella cinetica e di quella elastica
e dunque:
=1
2
Z
0
((
)2 +(
2
2)2)
Inserendo allora lo sviluppo (104) segue:
=1
2
∞X=1
0
0
Z
0
()
()
+1
2
∞X=1
Z
0
2()
2
2()
2
(108)
Ora integrando due volte per parti e sfruttando il fatto che:
4()
4 = ()
= 1 2
si ha: Z
0
2()
2
2()
2 =
Z
0
4()
4()
=
Z
0
()
()
Di conseguenza la (108) diventa:
=1
2
∞X=1
(0
0 +)(
()
()
)
In base alla proprietà di ortogonalità dei modi di oscillazione e alla =
2() la precedente si semplifica ulteriormente nella:
=1
2
∞X=1
((0())
2 + 2(())2)(
()
()
) (109)
che rappresenta l’energia istantanea che la trave possiede all’istante . Per
≥ , dove è un tempo oltre il quale cessa il carico che sollecita la trave, il
termine ((0())
2 + 2(())2) (vedi Appendice sulla trasformata di Fourier)
37
uguaglia il modulo al quadrato della trasformata di Fourier del secondo membro
della (107) calcolata in corrispondenza di = , ossia¯()
¯2. Posto allora:
() +Z +∞
−∞()
− = f()(()
()
)
f() +Z +∞
−∞f()
−
inserendo nella (109) segue:
=1
2
∞X=1
¯f()
¯2(
()
()
)
Cosicchè l’energia che il carico trasferisce al −modo è proporzionale a¯f()
¯2(
()
()
). Tale espressione, come deve essere, è invariante rispetto
alla trasformazione ()
→ ()
in quanto sia il numeratore che il denominatore
dipendono quadraticamente da ()
. L’indeterminazione sulla costante molti-
plicativa dei modi di oscillazione, non ha dunque nessun effetto sul termine¯f()
¯2(
()
()
).
3.3.1 Metodo di Rayleigh-Ritz per la trave
Nel caso in cui la trave non abbia caratteristiche omogenee lungo il suo asse,
nell’equazione (95) il coefficiente sarà una funzione di : = (). L’equazione
delle oscillazioni non forzate sarà quindi:
22 + 2()44 = 0 (110)
Cercando allora delle soluzioni del tipo (modi di oscillazione):
( ) = ()
si ottiene l’equazione agli autovalori:
44 =
2
2() (111)
Posto allora: 2 + 1
R 02(), + 2
2,() + 2
2() 0 la precedente
diventa:
44 = () (112)
che, assieme alle condizioni al contorno che definiscono i vincoli, definisce
il problema agli autovalori per la trave non omogenea. La presenza di una
() 6= 1 rende la soluzione analitica complessa. Si può allora applicare il
metodo di Rayleigh-Ritz, cercando gli estremali del seguente funzionale:
[] +R 000()
00()R
0()2()
(113)
38
che può essere ottenuto moltiplicando entrambi i membri della (112) per
e integrando (va fatta una doppia integrazione per parti sul numeratore).
Le proprietà del funzionale (113) sono analoghe a quello già esaminato per i
problemi ellittici.
3.4 Equazione della piastra
Con piastra si intende un continuo in cui una dimensione, lo spessore, è molto più
piccola delle altre due dimensioni, ma che, a differenza della membrana, resiste
sia a taglio che a flessione. L’equazione che descrive le piccole deformazioni
( ) della piastra ortogonali al piano cui appartiene in condizioni di riposo
(piano ) è data dall’equazione del quarto ordine:
2
2+ 2∇4 = ( ) (114)
∇4(·) + ∇2∇2(·) = 4(·)4
+ 24(·)22
+4(·)4
dove:
=spostamento ortogonale al piano di appartenenza della piastra a riposo
( ) =forza per unità di area applicata all’areola centrata in ( )
=densità areolare=massa/area
2 = 3(12(1− 2))
=modulo elastico
=spessore
=coefficiente di Poisson
3.4.1 Deformazione della piastra sotto l’azione di un carico ripartito
Cominciamo col risolvere un problema statico, ossia quando la funzione e
quindi sono indipendenti dal tempo. L’equazione (114) diventa allora:
∇4 = + 2 (115)
Risolvere tale equazione vuol dire trovare la deformata della piastra sotto
l’azione del carico ripartito Supponiamo che la piastra a riposo occupi un
dominio piano Ω del piano e che le condizioni di vincolo siano condizioni di
incastro. Si tratta allora di risolvere il seguente problema:½ ∇4 = in Ω
= 0 = 0 su Ω(116)
Vediamo ora di dare una formulazione debole del problema (116). Sia allora:
ªValutiamo ora le funzioni che annullano la variazione prima di Si ha:
= −
2= 0⇒ =
= [] (124)
Si ricava poi immediatamente che:
= 2
ZΩ
∇2∇2Ω
= 2
ZΩ
Ω
cosicchè le (124) diventano:ZΩ
∇2∇2Ω = []
ZΩ
Ω ∀ ∈M (125)
Ammettendo che ∈ 4(Ω), applicando due volte le formule di Green-Gauss
al primo membro della (125) e tenendo conto che ∈ A ∈M, segue:ZΩ
(∇4− [])Ω = 0 ∀ ∈M (126)
e quindi l’equazione di Eulero associata al funzionale (123) è:
∇4 = [] (127)
Tale equazione ci dice che gli estremali del funzionale (123) sono le autofun-
zioni del problema agli autovalori biarmonico (122) e, in corrispondenza di tali
funzioni, [] è proprio l’autovalore come mostra la (127). Poichè si può di-
mostrare che gli autovalori del problema biarmonico (122) sono numerabili e che
c’è un autovalore più piccolo, si tratta di trovare gli estremali del funzionale
(123) e valutare il più piccolo [] = 1 2 .
Dal punto di vista numerico, costruita l’approssimante:
=
X=1
∈M (128)
si tratta di trovare gli estremali della funzione:
(1 2 ) + [
X=1
] =(1 2 )
(1 2 )
(1 2 ) +X
=1
= ( )
(1 2 ) +X
=1
= ( ) (129)
42
con:
+ZΩ
∇2∇2Ω
+ZΩ
Ω
+ [1 2 ]
Poichè [] e quindi sono funzioni omogenee di grado zero, possiamo
minimizzare col vincolo che = 1, ossia imporre che:
( + e(1−)) = − e = 0 = 1 2
da cui segue il problema agli autovalori:
= e = 1 (130)
Calcolato il più piccolo autovalore e1, questo approssimerà il vero 1 del
problema (122). Dal corrispondente autovettore 1, si costruisce immediata-
mente il primo modo di oscillazione dall’approssimante (128).
Utilizzando come funzioni () () per una piastra rettangolare di lati
, i modi di oscillazione di una trave incastrata ai due estremi di lunghezza
pari a per le () e per le (), con la solita notazione ad un indice,
per il caso di una piastra di lati = 1, = 15 con = 102 = 100,
si trova come primo autovalore approssimato: e1 = 72952 a cui corrisponde
l’autofunzione (primo modo di oscillazione) mostrato in Fig.15.
Fig.15: primo modo di oscillazione
In Fig.16 e 16bis sono riportati rispettivamente i modi di oscillazione: sec-
ondo, terzo, quarto e quinto corrispondenti rispettivamente ai seguenti autoval-
ori approssimati:
(174019 437489 442918 637597)
43
Fig.16: secondo e terzo modo di oscillazione
Fig.16bis: quarto e quinto modo di oscillazione
44
4 CENNI SUL METODO AGLI ELEMENTI
FINITI (FEM)
4.1 Considerazioni generali
Il metodo agli elementi finiti si basa, dal punto di vista teorico, sulla formu-
lazione debole di un dato problema differenziale, mentre, dal punto di vista com-
putazionale, su una duplice approssimazione. La prima è relativa al dominio
Ω dove è definito il problema differenziale. Tale dominio viene approssimato
nella sua geometria da un ricoprimento Ω fatto da tanti elementi finiti (trian-
goli ad es. ma anche altre forme, con =massimo diametro dei vari elementi).
All’aumentare del numero degli elementi e alla corrispondente diminuzione delle
loro dimensioni ( → 0), il dominio Ω viene sempre meglio approssimato dal
dominio Ω definito dall’unione dei vari elementi finiti. Il dominio Ω viene an-
che denominato mesh. La seconda approssimazione riguarda la soluzione della
formulazione debole del problema, che viene ottenuta mediante la costruzione di
un approssimante appartenente ad un insieme ad un numero finito di dimensioni
ed applicando infine il metodo di Galerkin.
4.2 Costruzione generale della procedura numerica
Supponiamo di avere la seguente formulazione debole di un dato problema:
[ ] = []∀ ∈M (131)
consistente nel trovare ∈ A con [ ] forma bilineare simmetrica definitapositiva (vedi (31)). La (131) va risolta nei seguenti spazi infinito dimensionali:
vettoriale M . La matrice + [ ] associata a tali piecewise è ovvi-
amente una matrice sparsa (molti zeri con qualche "raro" elemento diverso da
zero) in quanto essendo[ ] un integrale associato al prodotto di due piece-
wise (o meglio alle loro derivate) è diverso da zero solo se i nodi e sono
"sufficientemente vicini" in modo da non rendere identicamente nullo il prodotto
delle piecewise.
Due funzioni piecewise associate ad es. alle zone colorate nella figura seguente
danno ovviamente contributo nullo.
Lavorando con matrici sparse, si devono adottare opportuni algoritmi che
ottimizzino i tempi di calcolo riarrangiando opportunamente le matrici per "rag-
gruppare" gli elementi diversi da zero e lavorare con meno numeri.
Infine, senza entrare troppo nei dettagli, per rappresentare la funzione 0()
che verifica le condizioni alla Dirichlet nella (133) si costruiscono delle piecewise
lineari associate ai nodi della frontiera ⊆ Ω, piecewise che assumono
su tali nodi i valori prescritti dalla condizione di Dirichlet (funzione ). Sui
nodi (se ci sono) appartenenti a Ω − in cui sono assegnate eventuali
condizioni alla Neumann, essendo queste "inglobate" nella formulazione debole
(condizioni naturali), le piecewise assumono su tali nodi valori indeterminati
()
( = 1 2 ) che si aggiungono alle incognite ( = 1 2 ).
In tal caso ci sono ovviamente altre piecewise di frontiera a coefficienti inde-
terminati che aumentano la dimensione diM da a + . La 0() viene
quindi rappresentata come combinazione lineare di tali piecewise di frontiera,
con coefficienti della combinazione non tutti noti se vi sono anche condizioni
alla Neumann. Mathematica 10.2 supporta con estrema facilità il metodo FEM
per equazioni differenziali di tipo lineare.
49
5 METODI ALLE DIFFERENZE FINITE
5.1 Considerazioni generali
In quel che segue faremo riferimento ai soli metodi numerici basati su una ap-
prossimazione alle differenze finite.
Per schema numerico intendiamo pertanto una opportuna trasformazione di
una o più equazioni differenziali in un insieme di equazioni algebriche che coin-
volgono i valori che la o le funzioni incognite assumono in punti computazionali
detti nodi, definiti da un reticolo di discretizzazione che ricopre il dominio di
definizione delle variabili indipendenti.
Le considerazioni che seguono, si riferiscono soprattutto alle applicazioni di
tali schemi alle equazioni differenziali alle derivate parziali del primo ordine di
tipo iperbolico e anche parabolico. L’equazione tipo su cui lavoreremo avrà
dunque la forma:U
+
U
= E (136)
oppure se ammette forma conservativa la:
U
+
F(U)
= E (137)
dove:
U + [1( ) 2( ) ( )]
E + [1( U) 2( U) ( U)]
F(U)
U=
ed essendo una matrice × reale che, nel caso di sistemi iperbolici, ha
autovalori tutti reali ed è diagonalizzabile.
Il caso più semplice della equazione (136) è la cosiddetta scalar linear wave
equation che ha la forma:
+
= 0 (138)
dove ( ) è l’incognita e è una costante.
Tale equazione verrà spesso usata per mettere in luce in maniera semplice ed
intuitiva alcune proprietà degli schemi numerici. L’utilizzo e l’implementazione
di specifici codici, non può prescindere da una analisi delle proprietà che uno
schema numerico deve possedere per essere in grado di rappresentare in maniera
adeguata una forma approssimata di una assegnata equazione differenziale. Tre
sono le caratteristiche fondamentali che deve possedere uno schema numerico:
esso deve risultare consistente, convergente e stabile. Nel diagramma di
seguito riportato, si evidenzia come la consistenza definisca una relazione fra la
struttura dell’equazione differenziale e la struttura della corrispondente formu-
lazione discreta, la stabilità una relazione tra la soluzione numerica computata
50
e l’esatta soluzione dell’equazione discretizzata ed infine la convergenza una re-
lazione fra la soluzione esatta dell’equazione discretizzata e la soluzione esatta
dell’equazione differenziale.
⇓
struttura dell’equazione differenziale
mstruttura della corrispondente formulazione discreta
⇓
soluzione numerica computata
msoluzione esatta dell’equazione discretizzata
⇓
soluzione esatta dell’equazione discretizzata
msoluzione esatta dell’equazione differenziale
5.2 Analisi degli schemi numerici
5.2.1 Consistenza
Cominciamo con l’introdurre il concetto di consistenza con un semplice esempio.
Consideriamo l’equazione (138) che qui riscriviamo nella forma:
+ = 0 = cost 0 (139)
Posto allora:
+ (∆ ∆)
con ∆∆ 0 e interi, approssimando le derivate della nel generico
punto = ∆ = ∆ nel seguente modo:
' +1 −
∆ '
−
−1∆
(140)
51
e sostituendo nella (139) si ottiene:
+1 −
∆+
−
−1∆
= 0 (141)
che rappresenta l’equazione alle differenze associata alla (139) e alle formule
(140) adottate per approssimare le derivate. La (141) lega i valori che la
assume nei tre nodi N1,N2,N3 mostrati in Fig.2:
xj )1( xj
tn
tn )1(
x
t
N1 N2
N3
Fig.2
Sfruttando lo sviluppo in serie di Taylor di punto iniziale = ∆ = ∆
possiamo scrivere:½+1 =
+ | ∆+ | ∆22! + | ∆33! +
−1 =
− | ∆+ | ∆22!− | ∆33! + (142)
che sostituite nell’equazione alle differenze (141) fornisce:
(| + | )+h | ∆2!− | ∆2! + | ∆23! + | ∆23! +
i= 0
(143)
Tale relazione mette in luce un aspetto molto importante che lega l’equazione
alle differenze (141) e l’equazione differenziale (139). Si nota infatti che i termini
entro le parentesi tonde esprimono esattamente ciò che l’equazione differenziale
(139) impone nel punto ( ), mentre i termini entro le parentesi quadre
rappresentano dei termini "spuri" che rappresentano l’errore strutturale che si
commette nel passare dalla forma differenziale alla forma discreta. In generale,
52
tale termine prende il nome di errore di troncamento relativo al passaggio da una
equazione differenziale(U) = 0 ad una forma discretizzata(U ) = 0.
Alla luce di quanto detto, possiamo dare la seguente definizione.
Data una forma approssimata alle differenze finite (U ) = 0 di una
equazione differenziale (U) = 0, definiamo errore di troncamento nel
generico nodo ( ) il valore:
+¯(U
)− (U)|= =¯
Uno schema numerico si dice allora consistente se l’errore di troncamento
tende a zero al tendere a zero degli intervalli di discretizzazione ∆∆. In
simboli:⇔ lim
(∆∆)−0 = 0
Come si vede dalla (143) lo schema numerico definito dalla (141) è consis-
tente, essendo = (∆) + (∆). Sulla base di quest’ultima, lo schema
viene inoltre detto del prim’ordine nello spazio e nel tempo (se fosse stato
= (∆2) + (∆) avremmo detto del secondo ordine nello spazio e
del primo nel tempo).
Il concetto di consistenza è molto importante in quanto costituisce una prima
condizione necessaria affinché uno schema numerico possa fornire soluzioni che
approssimino bene quanto si vuole la soluzione esatta di una equazione differen-
ziale.
5.2.2 Convergenza
Discutiamo ora più in dettaglio il concetto di convergenza.
In termini molto intuitivi, diremo che lo schema numerico è convergente
alla soluzione esatta dell’equazione differenziale che discretizza, se la soluzione
esatta dell’equazione discretizzata approssima sempre meglio la soluzione esatta
dell’equazione differenziale al tendere a zero delle dimensioni del reticolo di
discretizzazione.
Formaliziamo quanto detto. Sia (U) = 0 l’equazione differenziale da
risolvere con assegnate condizioni al contorno e iniziali per = 0. Ad un certo
istante 0 siaU( ) la soluzione esatta della(U) = 0 in eU quella
esatta della (U ) = 0 nel punto = ∆ = ∆ ottenuta con identiche
condizioni iniziali e al contorno dell’equazione differenziale. Se si verifica che:
Convergenza⇔ lim(∆∆)−0
max
¯U(∆ ∆)−U
¯= 0 ,∀ 0
diremo che lo schema numerico è convergente. In generale non è necessario
per la convergenza che (∆∆) tendano a zero in maniera arbitraria in quanto
in genere il loro rapporto deve rispettare la condizione di stabilità più avanti
definita. E’ evidente inoltre che al tendere di ∆ a zero, il numero di passi
per raggiungere un fissato istante , tende all’infinito.
53
5.2.3 Stabilità
Premettiamo che l’analisi di stabilità che tratteremo, è una analisi lineare, nel
senso che a rigore si applica ad equazioni di tipo lineare. L’analisi di stabilità
di tipo non lineare, è a tutt’oggi non ancora sufficientemente sviluppata per
poter essere sintetizzata in termini sufficientemente generali. Va però detto che
l’analisi lineare fornisce informazioni estremamente utili se utilizzate come linee
guida per i problemi non lineari, ferma restando la necessità di confermarne la
validità in campo non lineare mediante sperimentazione numerica.
Consideriamo allora l’equazione (136) supponendo che la matrice e il vet-
tore E siano costanti, in modo da considerare una equazione lineare. Sia poi
U( ) una soluzione della equazione (136) eW( ) una perturbazione sovraim-
posta a U. Imponendo allora che U( ) +W( ) soddisfi ancora la (136) e
ricordando che U( ) per ipotesi la soddisfa, segue cheW( ) deve verificare
l’equazione:W
+
W
= 0 (144)
ossia l’omogenea associata alla (136). Studiando come evolve la pertur-
bazione W, si analizza la stabilità della soluzione U( ) rispetto ad una per-
turbazione sovraimposta. Si può dimostrare, data la linearità della (144), che
per lo studio della stabilità è sufficiente analizzare il comportamento di soluzioni
elementari della (144) del tipo:
W = V(−) (145)
dove +√−1 un fissato numero d’onda reale e V uno scalare e un
vettore da determinare. Imponendo che la W soddisfi la (144), si ricava im-
mediatamente che (V) devono essere rispettivamente una coppia autovalore-
autovettore della matrice Essendo per ipotesi il sistema (144) iperbolico,
segue che gli autovalori di sono tutti reali e quindi anche lo è. Il significato
fisico di è chiaramente quello di una celerità dell’onda sinusoidale rappresen-
tata dallaW . Segue che il modulo della perturbazione elementareW risulta
costante al variare di e pari a:
|W | = |V|¯(−)
¯= |V|
in quanto se è reale¯(−)
¯= 1. Poichè la perturbazione W non si
amplifica nel tempo, possiamo dire che la U( ) è una soluzione stabile.
Sia ora dato uno schema numerico (U ) = 0 che discretizzi la (136),
e sia (U ) = 0 lo schema numerico che si ottiene ponendo E→ 0 in
(U ), ossia lo schema numerico che discretizza l’omogenea associata (144).
L’analisi di stabilità di tale schema viene fatta paragonando il comportamento
della soluzione elementare (145) con il comportamento di soluzioni elementari
analoghe ma relative all’equazione discretizzata (U ) = 0. Per la mag-
gior parte degli schemi di interesse applicativo (diremo più avanti quali sono
questi schemi), tali soluzioni elementari possono essere poste nella forma:
W = V(∆−∆) (146)
54
dove V è lo stesso vettore che compare nella (145) e ”” prende il nome di
celerità numerica. Anticipiamo che, a differenza di che è reale, in genere
è un numero complesso che si può quindi porre nella forma = + . Di
conseguenza segue:
W = V(∆−∆)∆ =⇒¯W
¯= |V| ∆ (147)
Di conseguenza se 0 la perturbazione cresce esponenzialmente nel tempo
e lo schema numerico è instabile, se ≤ 0 la perturbazione non cresce in moduloe lo schema è stabile. In pratica per valutare la stabilità di uno schema numerico
si procede nel seguente modo. Ponendo Y + V(−∆) la (146) diventa:
W = Y(∆) (148)
Imponendo che questa soddisfi la (W ) = 0, si arriva sempre ad
una relazione del tipo:
Y+1 = Y (149)
dove è una matrice a valori complessi detta di amplificazione che dipende
dallo schema numerico, da ∆ e ∆. Gli schemi numerici per i quali il vettore
V è lo stesso sia per la (145) che per la (146) (che sono peraltro quelli di maggior
interesse), sono quelli per i quali la matrice è un polinomio della matrice
ossia = 0 + 1+ 22 + , (con i coefficienti dipendenti da ∆ e ∆)
e di conseguenza V è autovettore sia di che di anche se in generale con
autovalori diversi.
L’azione lineare della matrice su Y rappresenta la trasformazione che
Y subisce nel passaggio discreto dall’istante ∆ all’istante (+1)∆. Poichè
dalla (148) segue che¯W
¯= |Y| porremo la seguente definizione: lo schema
numerico (U ) = 0 è stabile se:¯Y+1
¯≤ |Y|⇐⇒ |Y| ≤ |Y|⇒ Stabile (150)
viceversa: ¯Y+1
¯ |Y|⇐⇒ |Y| |Y|⇒ Instabile (151)
Lo studio della stabilità è dunque ricondotto alla determinazione delle con-
dizioni che devono soddisfare ∆ e ∆ affinchè sia soddisfatta la disequazione:
|Y| ≤ |Y| ∀Y ∈ C ∀ ∈ (152)
Si possono allora verificare tre casi:
1) la (152) è verificata ∀∆,∆ 0. In tal caso lo schema numerico si dice
linearmente incondizionatamente stabile.
2) la (152) non è verificata per nessun valore ∆,∆ 0. In tal caso lo
schema numerico si dice linearmente incondizionatamente instabile.
3) la (152) è verificata per quei valori di ∆,∆ 0 che verificano una
disuguaglianza del tipo (∆,∆) ≤ 0. In tal caso lo schema numerico si dice
55
linearmente condizionatamente stabile e la (∆,∆) ≤ 0 prende il nome di
condizione di stabilità.
L’importanza della stabilità di uno schema numerico è legata al cosiddetto
teorema di Lax. Poichè la consistenza da sola è una condizione unicamente
necessaria per la convergenza, il teorema di Lax asserisce che, almeno per le
equazioni lineari con condizioni iniziali e al contorno ben poste, consistenza e
stabilità di uno schema numerico, insieme, sono condizioni necessarie e suf-
ficienti per la convergenza della soluzione discretizzata alla soluzione esatta
dell’equazione differenziale.
5.3 Schemi espliciti e schemi impliciti
Gli schemi numerici alle differenze finite possono dividersi in due grandi cate-
gorie: schemi espliciti e schemi impliciti. Sia al solito (U ) = 0 una forma
discretizzata dell’equazione (U) = 0. Diremo che lo schema adottato è es-
plicito se può essere posto nella forma:
U+1 = F(UnUn−1 Un−k) (153)
intero ≥ 0
Il secondo membro della precedente indica simbolicamente che la dipendenza
funzionale è relativa ai valori che U assume negli istanti ∆ (−1)∆ (−)∆ calcolati in opportuni punti (−)∆ ()∆ (+)∆ con ≥ 0interi. Uno schema di questo tipo è detto a ( + 2)−livelli in quanto coinvolgegrandezze relative a ( + 2) istanti.
Sono impliciti gli schemi che non possono essere posti nella forma (153).Gli
schemi impliciti coinvolgono dunque almeno due valori che laU assume all’ultimo
livello ( + 1)∆. Questo fatto rende accoppiate le equazioni nelle incognite
U+1 , la cui determinazione richiede sempre la soluzione di un sistema di
equazioni algebriche. Negli schemi espliciti, invece, le equazioni possono essere
risolte una alla volta, essendo disaccoppiate. Un’altra importante differenza fra
schemi espliciti ed impliciti è che questi ultimi sono generalmente incondizion-
atamente stabili, mentre gli espliciti sono condizionatamente stabili e possono
essere utilizzati unicamente nel rispetto della condizione di stabilità che impone
determinate limitazioni al ∆ di integrazione.
5.4 Condizione di stabilità di Von Neumann
Un metodo molto usato per analizzare la stabilità di uno schema numerico,
è quello proposto da Von Neumann. Questo metodo si basa sulla analisi degli
autovalori della matrice di amplificazione . Siano 1 2 gli autovalori
(in genere complessi) di . Si definisce raggio spettrale della matrice (che
indicheremo con ()) il più grande fra i moduli degli autovalori: () =
max(|1| |2| | |). La condizione di stabilità di Von Neumann afferma checondizione necessaria per la stabilità di uno schema numerico è che () ≤ 1.Tale condizione diventa anche sufficiente se è diagonalizzabile. Per quegli
56
schemi numerici per i quali V è autovettore sia della matrice che di , la
diagonalizzabilità di implica quella , ed in tal caso la condizione di stabilità
di Von Neumann è necessaria e sufficiente.
Facciamo ora una semplice applicazione di tale metodo alla scalar wave equa-
tion. Tale applicazione permette inoltre di mettere in luce una interessante
interpretazione euristica della instabilità numerica.
Consideriamo allora l’equazione:
+ = 0 = cost 0 (154)
e discretizziamola come già visto:
+1 −
∆+
−
−1∆
= 0⇐⇒ +1 −
+ ( −
−1) = 0 (155)
con + ∆∆. Posto allora = ( + ∆) imponendo che
questa soddisfi la (155), si trova +1 = essendo il fattore di amplifi-
cazione (in tal caso la matrice è 1x1 e si parla di fattore) definito da:
= 1 + (− − 1)
Segue calcolando il modulo dopo qualche semplificazione:
|| =p1 + 2(1− cos())( − 1)⇒ || ≤ 1⇐⇒ ≤ 1
La condizione di stabilità per lo schema (155) è dunque:
∆∆ ≤ 1
In particolare per = 1 la (155) fornisce +1 =
−1. Tale relazione,tenendo conto che la soluzione esatta della (154) è una funzione del tipo =
(−), fornisce in tal caso una soluzione discreta che nei nodi computazionalicoincide con la soluzione esatta dell’equazione differenziale (154).
Diamo ora una spiegazione dell’instabilità numerica, sfruttando un’analogia
fisica. Nell’equazione discretizzata (155) inseriamo gli sviluppi di Taylor (142).
Sfruttando inoltre la relazione:
+ = 0⇒ = − = − = −(−) = 2
si arriva dopo alcuni passaggi alla:
(| + | ) =(1− )
2∆ | +(∆2 +∆2) (156)
L’equazione da cui siamo partiti "+ = 0" è una equazione iperbolica
di puro trasporto. Lo schema numerico (155) essendo equivalente alla (156)
impone (nel generico nodo computazionale), invece che una equazione di puro
trasporto iperbolica, una equazione di trasporto-diffusione di tipo parabolico a
causa della presenza del termine di derivata seconda in a secondo membro. Il
57
coefficiente(1−)2
∆ può essere interpretato, in analogia alle equazioni idrodi-
namiche, come una viscosità artificiale ”” che lo schema numerico introduce.
Rispetto però alla viscosità reale che è sempre positiva dovendo dissipare ener-
gia, la ”” è dotata di segno.
Se infatti 1 (condizione di stabilità) segue che 0 e il termine
diffusivo | diventa un termine"pozzo" che smorza l’onda durante il suopropagarsi. Se = 1 segue che = 0, non vi è dissipazione e lo schema
numerico fornisce la soluzione esatta. Se infine 1 (condizione di instabilità)
segue che 0. Una viscosità negativa rende il termine diffusivo un termine
"sorgente" ed in tal caso l’onda nel suo propagarsi si amplifica senza limite,
rendendo instabile lo schema numerico
Concludendo tale paragrafo, possiamo riassumere brevemente quanto ot-
tenuto: consistenza e stabilità di uno schema numerico, studiate in ambito
lineare, sono delle condizioni imprescindibili che ogni schema deve soddisfare
per poter essere applicato con successo. Come abbiamo visto, limitandoci a
equazioni lineari, tali condizioni sono necessarie e sufficienti per la convergenza
della soluzione approssimata fornita dallo schema numerico alla soluzione esatta
dell’equazione differenziale.
Nell’ambito delle equazioni non lineari, le suddette condizioni non sono in
generale sufficienti pur rimanendo necessarie. Allo stato attuale delle conoscenze
infatti, solo la sperimentazione numerica può "empiricamente confermarne" la
sufficienza.
5.5 Dissipazione e Dispersione
Ritorniamo al confronto fra le soluzioni elementari della (144) e quelle di un
associato schema numerico:
W = V(−)
W = V(∆−∆) = Y∆ (157)
Y + V(−∆)
Come abbiamo visto la W soddisfa l’equazione discretizzata se e solo se
Y+1 = Y. E quindi in base alla definizione di + V(−∆) se e solose (I matrice identità):
(− I−∆)V = 0 (158)
Dunque la V(∆−∆) soddisfa l’equazione discretizzata se e solo se Vè autovettore della matrice di amplificazione e −∆ è il corrispondenteautovalore. Assumendo quale incognita la celerità numerica complessa , per
ricavarla si deve risolvere l’equazione(−I) = 0 e ricavare dalla relazione: = −∆
Ovviamente avremo tante quanti sono gli autovalori di . Considerando
allora il − autovalore ( = 1 2 ) segue:
= −∆ ⇒ = − arg()(∆) + log(||)(∆) (159)
58
avendo sfruttato la definizione di logaritmo di un numero complesso per la
quale: log() + log(||)+ arg() dove arg(·) indica l’argomento del numerocomplesso, ossia l’angolo (positivo in verso antiorario) che il vettore che lo rap-
presenta forma nel piano complesso con l’asse reale. Sostituendo nella (157) e
ragionando sui − autovalori e autovettori segue:
W = V(−)
W = VIm()∆(∆−Re()∆) (160)
Re() = − arg()(∆) Im() = log(||)(∆)
Confrontando laW e laW , notiamo allora delle importanti differenze.
Innanzi tutto, il fatto che sia complessa implica, come era già stato messo
in evidenza, che il modulo della soluzione numerica non è costante ma variabile
nel tempo secondo un fattore esponenziale Im()∆. Quindi se Im()
0⇔ || 1 (instabilità) la soluzione numerica si amplificherà indefinitamenteal crescere di . Se invece Im() ≤ 0 ⇔ || ≤ 1 (stabilità) la soluzione
numerica non crescerà al crescere di , tendendo a smorzarsi nel caso valga la
disuguaglianza stretta e a rimanere costante in modulo nel caso di uguaglianza.
Diremo allora che uno schema numerico è dissipativo relativamente al −
autovalore di se || 1. Il livello di dissipazione di uno schema numericopuò essere controllato (reso arbitrariamente vicino all’unità || . 1) agendo
sull’ampiezza di ∆∆.
Un’ altra importante differenza che emerge dal confronto delle W e la
W , è che Re() è l’analogo di e rappresenta quindi la celerità con cui lo
schema numerico fa evolvere l’onda W . Va però tenuto conto che, essendo
per ipotesi in ambito lineare, è una costante (in particolare non dipende da
) mentre Re() = (∆∆ ). Quindi mentre le varie armoniche W ,
al variare di , viaggiano tutte alla stessa velocità , le armoniche W , al
variare di , viaggiano in generale a velocità differenti. In analogia con le onde
della fisica, tale aspetto della soluzione numerica prende il nome di dispersività
e lo schema numerico si dice in tal caso dispersivo. Anche la dispersività può
essere controllata agendo sull’ampiezza di ∆∆.
In letteratura la bontà di uno schema numerico viene spesso analizzata
facendo riferimento ai seguenti tre parametri:
+Re()
+ ∆( −Re()) + Im()∆)
detti rispettivamente rapporto dispersivo,errore di fase,coefficiente di ampli-
ficazione. Inoltre è abbastanza evidente che la consistenza implica → 1 →0 → 1, facendo tendere opportunamente a zero ∆ e ∆.
59
6 APPENDICE: Approfondimenti e Richiami
6.1 Formule di Green-Gauss
Sia Ω ⊂ un insieme aperto, n = (1 2 ) la normale esterna alla sua
frontiera Ω e (1 2 ) una funzione ∈ 1(Ω ∪ Ω). Allora valgono le:ZΩ
Ω =
ZΩ
Applicando le precedenti ad una funzione = prodotto di due funzioni si
ricava immediatamente la formula di integrazione per parti:ZΩ
Ω =
ZΩ
−ZΩ
Ω
6.2 Ortogonalità degli autovettori nei problemi generaliz-
zati agli autovalori
Consideriamo il problema generalizzato agli autovalori:
= (161)
dove e sono matrici simmetriche ad elementi reali di dimensione (×)e + (1 2 ) . Dimostriamo che detti12 due autovettori corrispon-
denti rispettivamente ad autovalori distinti 1 2, si ha:
1 2 = 0 (162)
La precedente, nel caso che sia la matrice identità, ci dice che autovettori
corrispondenti ad autovalori distinti, sono tra loro ortogonali.
Dim.
Per ipotesi: ½1 = 11
2 = 22(163)
segue:
(1) = 1(1)
⇒ 1 = 1
1 (164)
essendo per ipotesi = = . E quindi, rispettivamente per la(164)
e per la (163), si ha:½1 2 = (
1 )2 = 1
1 2
1 2 =
1 (2) = 21 2
⇒ 11 2 = 2
1 2
da cui:
(1 − 2)1 2 = 0
Essendo per ipotesi 1 − 2 6= 0⇒ 1 2 = 0 ossia la tesi.
60
6.3 Ricerca del minimo assoluto di come minimo vinco-
lato
Verifichiamo che se il funzionale [] = N []D[] definito dalla (42) gode dellaproprietà: [] = N []D[] = 2N []
2D[] = [], ∀ 6= 0 (ossia il funzionale èomogeneo di grado zero), tale proprietà garantisce che la ricerca del minimo di
[] in A0 − 0 è equivalente alla ricerca del minimo del numeratore N [] colvincolo che il denominatore D[] valga uno.Dim
Sia 1 ∈ A0 − 0 la funzione sulla quale [] assume il minimo assoluto inA0 − 0. Dalla [] = [] segue che ∀ 6= 0 anche su 1 ∈ A0 − 0, []assume il minimo assoluto. Ma [1] = N [1]D[1] = 2N [1]
2D[1] . Scegliendo
allora = 1pD[1] segue che:
[1] = N [1]1e dunque il minimo di [] viene assunto su una funzione 1 che rende
unitario il denominatore. Da ciò segue che:
min∈A0−0
[] = min∈A0D[]=1
N []
6.4 Ricerca degli autovalori successivi al primo col metodo
di Rayleigh-Ritz
Ricordando la (44), segue che gli autovettori e gli autovalori approssimati succes-
sivi al primo, diciamo il − autovettore e autovalore con = 2 3 possono
essere trovati minimizzando = (vedi (54)) col vincolo 1 − =
1− = 0 e con gli ulteriori vincoli: = 0, = 1 2 − 1, dove
sono gli autovettori "precedenti" a . Utilizzando il metodo dei moltipli-
catori di Lagrange segue che bisogna risolvere le + equazioni:⎧⎪⎪⎨⎪⎪⎩ = 0, = 1 2
= 1
= 0 = 1 2 − 1
+ + e(1−) +P−1
=1
(165)
nelle + incognite + (1 2 ) ,1,2,,−1,e, (essendo questiultimi i moltiplicatori di Lagrange). Sia ora (e) la coppia −
autovalore-autovettore (con e1 e2 e) del problema = e verifi-
cante la = 1. Verifichiamo allora che le + grandezze ( = 0,
= 1 2 − 1), e verificano le (165). Inserendo infatti in queste equazionitali grandezze segue:
= ( + e(1−)) = 0⇒ = e
= 1
= 0 = 1 2 − 1
61
Le prime due equazioni sono verificate per definizione di (e). Le restanti
− 1 equazioni sono verificate in base alla proprietà di ortogonalità (162). Inconclusione, la ricerca dei minimi vincolati di per l’approssimazione degli
autovalori successivi al primo, in virtù della proprietà di ortogonalità (162),
si ottiene semplicemente risolvendo il problema agli autovalori = ee ordinando in ordine crescente gli autovalori approssimati e1 e2 e .Allora ogni e ( = 1 2 ) approssimerà il valore vero e ogni permetterà
di costruire la corrispondente autofunzione sulla base dell’approssimante (53).
6.5 Trasformata di Fourier: interpretazione energetica
Consideriamo una funzione () con ∈ (−∞+∞) che verifichi la condizione:Z +∞
−∞| ()| +∞ (166)
Si definisce allora trasformata di Fourier della (), la funzione () così
definita ( +√−1):
() +Z +∞
−∞ ()− (167)
La () è una funzione a valori complessi della variabile reale . Tale
variabile può ovviamente essere interpretata come una pulsazione. Vediamo ora
quale interpretazione "energetica" può essere data alla (). Consideriamo
un oscillatore armonico forzato e non smorzato di massa unitaria descritto dalla
equazione:
00() +2() = () (168)
dove () (pensata come forza eccitatrice) è una funzione che supporremo
nulla per 0 e che verifica la condizione di integrabilità (166). Supponiamo
inoltre che all’istante = 0 l’oscillatore sia fermo (0(0) = 0) e in posizione diriposo ((0) = 0).
Avendo supposto l’oscillatore di massa unitaria, la sua energia totale istan-
tanea, somma della cinetica e potenziale, sarà:
=1
2(0())2 +
1
22(())2 (169)
Associamo ora alla funzione (), la funzione a valori complessi () così
definita:
() + 0() + () (170)
il cui modulo al quadrato moltiplicato per 12, fornisce l’energia dell’oscillatore:
=1
2|()|2 (171)
Segue, derivando rispetto al tempo ():
0() + 00()+
0() =
00()+(()−()) =
00()+2()+()
62
e quindi tenendo conto della (168), segue che la () verifica l’equazione
differenziale del prim’ordine:
0()− () = () (172)
Le condizioni iniziali (0) = 0 0(0) = 0 inducono la condizione iniziale
(0) = 0(0) + (0) = 0. Moltiplicando allora la (172) per − segue:
0()− − ()− = ()− ⇒
(()−) = ()−
Integrando allora fra = 0 e un generico istante 1 segue, tenendo conto che
(0) = 0:
(1)−1 =
Z 1
0
()−⇒ |(1)| =¯Z 1
−∞()−
¯avendo tenuto conto che la forza eccitatrice () (per ipotesi) è nulla per
0 e che essendo reale¯−1
¯= 1. Facendo tendere 1 → +∞ segue:
|(+∞)|2 =¯Z +∞
−∞()−
¯2=¯ ()
¯2= 2(+∞)
La precedente relazione suggerisce una interpretazione estremamente efficace
circa il significato della trasformata di Fourier () di una funzione () che
verifica la (166) ed è nulla per 0. Pensando infatti tale funzione come forza
eccitatrice di un oscillatore armonico di massa unitaria e pulsazione propria ,
che per = 0 è in posizione di riposo e fermo, il modulo al quadrato della trasfor-
mata di Fourier () della (), calcolato in corrispondenza della pulsazione
propria dell’oscillatore, fornisce il doppio dell’energia totale che la forza ()
trasferisce all’oscillatore nell’intervallo di tempo [0+∞).Consideriamo in particolare il caso realistico in cui la forza () è diversa
da zero in un intervallo di tempo finito. Questo significa che, forze eccita-
trici comunque diverse sia per durata che per forma, se in corrispondenza di
= hanno trasformate di Fourier con uguale modulo, avranno trasferito
all’oscillatore la stessa energia una volta terminata la loro azione.
6.6 Principio di minima azione
Il principio di minima azione di Hamilton, afferma che in un sistema mecca-
nico soggetto a forze derivabili da un potenziale − (e quindi dotato di energia
potenziale ), i moti effettivi del sistema sono gli estremali del seguente fun-
zionale:
=
Z 2
1
( − ) (173)
dove 1 e 2 sono due istanti di tempo entro i quali evolve il sistema ed
è l’energia cinetica. Il funzionale (173) prende il nome di azione, mentre il
63
principio si chiama di minima azione perché in molti problemi della meccanica,
l’estremale è effettivamente un minimo. Per poter derivare correttamente da
tale principio le equazioni della meccanica, è necessario specificare lo spazio cui
devono appartenere le variazioni degli argomenti del funzionale (173) (insieme
M) al fine di poter calcolare la variazione prima dell’azione. Tali variazioni
devono sempre annullarsi in corrispondenza di 1 e 2 I moti variati devono
quindi (nello spazio delle fasi del sistema) partire e arrivare negli stessi punti
del moto effettivo, come mostrato schematicamente in Fig.1 dove la linea più
spessa simboleggia il moto effettivo del sistema mentre le altre possibili moti
variati. Inoltre, qualora gli argomenti dell’azione dipendano anche da variabili
spaziali appartenenti ad un dominio Ω tali variazioni devono annullarsi anche
su Ω in modo che le funzioni variate verifichino le stesse condizioni al contorno
delle funzioni non variate. Nel simbolismo finora usato + è la funzione
variata, la variazione, la funzione non variata.
1t
2t
Fig.1
6.7 Deduzione dell’equazione della membrana dal princi-
pio di minima azione
Deduciamo ora dal Principio di Minima Azione, a titolo di esempio, l’equazione
della membrana (58). L’energia cinetica della membrana è data da:
=
ZΩ
1
2(
)2 =
ZΩ
1
2(
)2Ω (174)
L’energia potenziale associata alla deformazione elastica può essere scritta
come:
= 0(Area membrana deformata-Area membrana a riposo) ossia:
= 0(
ZΩ
s1 + (
)2 + (
)2Ω−
ZΩ
1Ω)
dove 0 è la tensione della membrana. Per piccole deformazioni si ha che:
64
s1 + (
)2 + (
)2 − 1 ' 1
2((
)2 + (
)2)
cosicchè risulta:
= 0(
ZΩ
1
2((
)2 + (
)2)Ω) (175)
Infine il potenziale associato ad un carico distribuito ( ) (carico in
direzione verticale) è:
= (176)
cosicchè il potenziale complessivo − = −+ . L’azione per il sistema
membrana è dunque:
[] =
Z 2
1
( − + ) =Z 2
1
ZΩ
∙1
2(
)2− 1
20((
)2 + (
)2) +
¸Ω (177)
Calcoliamo allora la variazione dell’azione. Per quanto detto sui moti variati:
