Notes 1819
Notes 1819
I Metodi numerici 1
I
2.4.2 Metodo di Newton come iterazione di punto fisso . . . . . . 44
2.5 Criteri di arresto . . . . . . . . . . . . . . . . . . . . . . . . . . . . 45
II
5.3.4 Theta-Metodi . . . . . . . . . . . . . . . . . . . . . . . . . . 114
5.3.5 Risoluzione di sitemi differenziali ordinari . . . . . . . . . . . 115
5.3.6 Risoluzione numerica di problemi del secondo ordine . . . . 116
6 Introduzione 121
6.1 Equazioni alle derivate parziali . . . . . . . . . . . . . . . . . . . . . 121
6.2 Condizioni al contorno . . . . . . . . . . . . . . . . . . . . . . . . . 122
6.3 Serie di Fourier . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 124
III
III Problemi in forma debole 185
IV
13.3 Problemi differenziali . . . . . . . . . . . . . . . . . . . . . . . . . . 257
13.3.1 Equazione di Laplace nel semipiano . . . . . . . . . . . . . . 257
13.3.2 Problema di diffusione . . . . . . . . . . . . . . . . . . . . . 259
V
Elenco delle figure
VI
14.3 Schema risoluzione ODE tramite trasformata di Laplace . . . . . . 269
VII
Elenco delle tabelle
VIII
Elenco dei codici MATLAB®
IX
X
Introduzione
Problema fisico Osservazione del fenomeno in natura, ipotesi sulle sue cause e
determinazione dei parametri che sono in grado di influenzarne il comporta-
mento.
Il risultato di questo lavoro è una soluzione affetta da più o meno grandi appros-
simazioni, che derivano in parte dalle semplificazioni apportate ai modelli mate-
matici e numerici e in parte dai processi calcolo. Al netto di queste considerazioni,
la soluzione trovata ha una grande validità in campo scientifico e ingegneristico;
per questo motivo negli ultimi anni le tecniche di calcolo numerico sono divenute
sempre più fondamentali e tuttora rappresentano molto spesso l’unica soluzione
disponibile alla risoluzione di problemi molto complessi.
XI
Problemi ben posti
Prima di iniziare ad analizzare le procedure numeriche è bene introdurre alcune
terminologie utili a classificare i metodi.
Definizione 1 (Problema ben posto). Un problema numerico si dice ben posto, se
esiste una sola soluzione che dipende con continuità dai dati.
Tale proprietà garantisce che a fronte di una piccola variazione dei dati in in-
gresso non si abbiano delle divergenze incontrollate della soluzione numerica, ma
che vari in modo continuativo. Per identificare in termini quantitativi la perturba-
bilità del sistema si fa uso di un parametro, chiamato numero di condizionamento.
Definizione 2 (Numero di condizionamento). Il numero di condizionamento è un
parametro in grado di stimare la stabilità di un metodo numerico ed è definito
come: { ∥δx∥ }
/∥x∥
k = sup ∥δd∥
δd /∥d∥
dove x è la soluzione e d il dato in ingresso.
Inoltre è utile definire una proprietà che metta in relazione i problemi mate-
matici con i problemi numerici.
Definizione 3 (Problema numerico consistente). Un problema numerico viene
detto consistente quando converge al problema matematico:
XII
Tabella 1: Principali tipologie di convergenza ordinate per velocità in modo crescente.
Definizione 5 (Insieme dei Floating Point F). Si chiama insieme dei Floating
Point F(β, t, L, U ) l’insieme dei numeri che possono essere scritti come:
x ∈ F ⇐⇒ x = (−1)s mβ e−t
XIII
Da questa definizione risoluta chiaro che la minima differenza percettibile dal
computer varia a seconda del numero, anche per questo è definito un parametro.
.
εm = min{1 + x ̸= 1}
x∈F
È importante notare che per definire un insieme di numeri floating point è sempre
necessario definire la mantissa β, il numero delle cifre della mantissa t e i due
estremi L e U (entrambi sempre positivi).
XIV
Tabella 2: Principali formati per la visualizzazione del numero π all’interno
dell’ambiente MATLAB® (* opzione di default).
format comando
XV
Parte I
Metodi numerici
1
Capitolo 1
Ax = b
∃! x0 ∈ R ⇐⇒ det A ̸= 0
x0 = A−1 b
Purtroppo il calcolo della matrice inversa A−1 , necessita della matrice dei cofattori,
che a sua volta si ottiene tramite il determinante di ognuna delle sottomatrici che
in generale è molto lungo e complesso.
Un’alternativa consiste nell’applicare il metodo di Cramer, che riesce a calcolare
2
le componenti del vettore x tramite il rapporto fra il determinate della matrice Bk
e il determinante della matrice A (dove la matrice Bk è ottenuta sostituendo la
k-colonna della matrice A con il vettore dei temrini noti b).
det(Bk ) [ ]
xk = Bk = a1 | . . . | ak−1 | b | ak+1 | . . . | an
det(A)
∑
n
det(A) = (−1)i+j aij det(Aij ) ∀i = 1, 2, . . . , n
j=1
Metodi diretti Metodi che arrivano alla soluzione esatta del sistema in un nu-
mero finito di passaggi.
Metodi iterativi Metodi che giungono alla soluzione esatta solamente dopo un
numero infinito di passaggi; ovviamente ci si può accontentare di una solu-
zione approssimata.
Ovviamente negli ambienti software che si usano tutti giorni queste tecni-
che sono già implementate e non serve ogni volta scriverle da capo: ad esempio
3
in MATLAB® è disponibile la funzione backslash che automaticamente sceglie il
metodo più efficiente per risolvere il sistema lineare:
x = A\b
Nello specifico A è la matrice, b il vettore dei termini noti e x la soluzione del si-
stema; essendo una moltiplicazione fra matrici è importantissimo che A moltiplichi
da sinistra.
Dx = b
∏
n
det(D) = dii ̸= 0 ⇐⇒ dii ̸= 0 ∀i = 1, 2, . . . , n (1.2)
i=1
4
1.1.2 Risoluzione di sistemi triangolari
La seconda tipologia di sistemi facili da risolvere è costituita dai sistemi triangolari
(inferiori o superiori); tali nomi derivano dalla forma che caratterizza la matrice
A, cioè quella di un triangolo.
Ossia quando gli elementi della matrice A sopra la diagonale principale sono nulli.
Nonostante siano presenti molti più termini non nulli rispetto alla matrici dia-
gonali (Paragrafo 1.1.1), è ancora possibile calcolare il determinante tramite la
moltiplicazione di tutti i termini presenti sulla diagonale e vale:
∏
n
det(A) = aii ̸= 0 ⇐⇒ aii ̸= 0 ∀i = 1, 2, . . . , n (1.3)
i=1
Per risolvere questi sistemi viene solitamente adottato l’algoritmo della sosti-
tuzione in avanti (o forward substitution), leggermente più complesso rispetto al
precedente.
L’idea è quella di risolvere la prima equazione, la cui soluzione è immediata per-
ché è disaccoppiata da tutte le altre. A questo punto si può sfruttare l’informazione
appena trovata per disaccoppiare e risolvere la seconda e così via. . . l’algoritmo
finale risulta quindi:
b1
x1 =
a11
Per i = 2, 3, . . . , n
( ∑
i−1 )
1
xi = bi − aij xj
aii j=1
Per risolvere il sistema sono quindi necessarie n divisioni (una per ogni riga)
ed un numero di somme e sottrazioni che è pari a zero al primo passaggio aumenta
5
di un’unità ad ogni incremento di i:
∑n
n(n − 1)
#flops = n + 2 (i − 1) = n + 2 = n2 (1.4)
i=1
2
Ossia quando gli elementi della matrice A sotto la diagonale principale sono nulli.
Ancora si disaccoppiano una a una le equazioni, ma in questo caso si inizia con
l’ultima equazione e si termina con la prima. L’algoritmo necessario per questa
operazione prende il nome di sostituzione all’indietro (o backward substitution) ed
è definito nel modo seguente:
bn
xn =
ann
Per i = n − 1, n − 2, . . . , 1
( ∑n )
1
xi = bi − aij xj
aii j=i+1
Dato che questo procedimento è del tutto analogo al precedente, si può facil-
mente dimostrare che il numero di operazioni necessarie al completamento della
risoluzione è O(n2 ) (crf. Eq. (1.4)).
6
si fa largo uso del metodo di fattorizzazione LU, che consiste nel costruire due ma-
trici L e U (con L, U ∈ Rn×n ) tali che A = LU. In particolare si ricercano le due
matrici facendo in modo che L sia triangolare inferiore e U triangolare superiore;
una volta trovate si possono risolvere i sistemi:
Ly = b per trovare y
Ax = L(Ux) = Ly = b −→ (1.5)
Ux = y per trovare x
Una definizione simile può essere scritta per le matrici definite negative sem-
plicemente sostituendo il maggiore stretto con un minore stretto all’interno della
definizione. Stesso discorso vale per le matrici definite semipositive o seminega-
tive, in cui si trova nella definizione una disuguaglianza semplice al posto di una
disuguaglianza stretta.
∑
n
|aii | ≥ |aij | ∀i = 1, 2, . . . , n
j=1
j̸=i
∑
n
|aii | ≥ |aji | ∀i = 1, 2, . . . , n
j=1
j̸=i
7
Teorema 1.1.1 (Condizione necessaria e sufficiente per fattorizzazione LU). Data
una matrice A ∈ Rn×n non-singolare esiste una sola fattorizzazione LU se e solo
se le sottomatrici principali di nord-ovest Ai sono tutte non-singolari.
∃! A = LU ⇐⇒ det(Ai ) ̸= 0 ∀i = 1, 2, . . . , n
Si può facilmente notare che si sta cercando di determinare sei elementi con sola-
mente quattro equazioni disponibili. In particolare eseguendo il prodotto matriciale
si ottiene il seguente sistema sotto determinato:
a11 = l11 u11
a = l u
12 11 12
a = l u
21 21 11
a = l u + l u
22 21 12 22 22
8
Questo ragionamento può essere esteso a una matrice generica A ∈ Rn×n : in questo
caso si hanno a disposizione n2 vincoli, per determinare n2 + n incognite. Si ha
quindi la possibilità di determinare a priori n elementi delle matrici L e U; per
convenzione si assume che:
lii = 1 ∀i = 1, 2, . . . , n
Si ottiene quindi una matrice L che non solo è diagonale inferiore, ma che ha anche
ha tutti gli elementi pari a 1 sulla diagonale principale. Sotto questa ipotesi si può
dimostrare che le matrici L e U vengono calcolare con il metodo di eliminazione di
Gauss (MEG.), in particolare, ad ogni iterazione si riarrangia la matrice A al fine
di farla diventare una matrice triangolare superiore, determinando volta per volta
una colonna di L.
in cui si può ben distinguere il blocco non ancora lavorato di termini non nulli
in basso a destra e la matrice diagonale superiore finita, che avanza da sinistra a
destra ad ogni iterazione.
(k)
Definizione 1.1.7 (Elementi pivotali). Tutti gli elementi akk della matrice A(k)
dell’Eq. (1.6) si chiamano elementi pivotali.
Se tutti gli elementi pivoltali della matrice generica A ∈ Rn×n sono non nulli,
9
allora si può definire l’algoritmo del MEG.:
A(1) = A
Per k = 1, 2, . . . , n − 1
(k)
aik
lik = (k)
∀i = k + 1, . . . , n
akk
(k+1) (k) (k) ∀i = k + 1, . . . , n
aij = aij − lik akj
∀j = k, . . . , n
È interessante notare che per calcolare gli elementi lik sono necessari due cicli an-
(k+1)
nindati, mentre per sistemare gli elementi aij ne servono tre. Si può dimostrare
che per eseguire questo tipo di lavoro sono necessari un numero di cicli dell’ordi-
ne di O(n3/3), di conseguenza le sostituzioni in avanti e all’indietro necessarie per
risolvere il sistema lineare diventano trascurabili. Infine, se si deve lavorare con
una matrice A è di dimensione molto elevata, è possibile, al fine di risparmiare
spazio occupato, sovrasciverla ad ogni passaggio, in modo che tutti gli elementi
sulla diagonale e sopra di essa siano gli elementi della matrice U, mentre quelli
al di sotto della diagonale siano i coefficienti della matrice L (questo è possibile
solamente perché tutti gli elementi diagonali della matrice L sono noti e pari a
uno).
Questo algoritmo è una variante del MEG. ed è utilizzato nel caso in cui, al fine di
determinare la fattorizzazione LU di una matrice, sia presente un termine pivotale
nullo. Si sa per certo che la matrice ha almeno un pivot diverso da zero ad ogni
passaggio, altrimenti la matrice stessa sarebbe singolare (evento che si è escluso
nelle ipotesi iniziali).
Il meccanismo alla base di questo metodo è quello di effettuare delle permu-
tazioni sulla matrice A, atte a scambiare le sue righe, in modo che sia sempre
possibile controllare l’elemento pivotale. Per fare ciò, basta moltiplicare A e il
vettore dei termini noti per una matrice P, detta di permutazione; a questo punto
è necessario rivedere la risoluzione dei sitemi presenti nell’Eq. 1.5, modificandoli
10
in questo modo:
Ly = Pb per trovare y
PAx = L(Ux) = Ly = Pb −→ (1.7)
Ux = y per trovare x
La tecnica del pivoting serve anche a ridurre la perdita di precisione nei calcoli.
Nello specifico si può dimostrare che conviene utilizzare l’elemento pivotale di
(k)
modulo massimo tra tutti quelli disponibili, amk :
(k) (k)
akk = max |amk |
m=k,...,n
A questo punto occorre scambiare la k-esima con la m-esima riga e per farlo ser-
va la matrice P adatta. Tale matrice è una matrice di zeri, su cui vengono inseriti
degli 1 nella posizione corrispondente di ogni elemento pivotale estratto. La ma-
trice P ha anche le proprietà di essere simmetrica, ortogonale e aver determinante
pari a uno.
P = PT = P−1 det(P) = 1
[L U P] = lu(A)
11
1.1.4 Algoritmo di Thomas
Questo algoritmo viene utilizzato per risolvere una particolare classe di sistemi
lineari: quelli caratterizzati da una matrice di tipo tridiagonale.
a1 c 1 0 ... 0
e 2 a2 c 2 0
0 e 3 a3 c 3 0
A=
.. .. ..
0 . . .
.
.. c
n−1
0 en an
Come si può facilmente notare, le matrici L e U sono molto più semplici rispetto
a prima, perché contengono molti più zeri, e per questa ragione si avranno gradi
semplificazioni sull’algoritmo:
α 1 = a1
Per i = 1, 2, . . . , n
ei
βi =
αi−1
αi = ai − βi ci−1
12
x = yn
n
Ux = y −→ αn
i = 1 [y − c x ] i = n − 1, n − 2, . . . , 1
i αi i i i+1
dove R è una matrice triangolare superiore, con elementi positivi sulla diagonale.
L’algoritmo necessario è quindi:
Per i = 1, 2, . . . , n
v
u
u ∑i−1
t
rii = aii − r2
ki
k=1
[ ∑
i−1 ]
1
rij = aij − rki rkj ∀j = i + 1, . . . , n
rii k=1
Il fatto di implementare una sola matrice fa risparmiare molta memoria (la metà
rispetto alla fattorizzazione LU). Purtroppo non si hanno grandi benefici per quan-
to concerne all’utilizzo delle risorse computazionali: infatti il costo del processo è
dell’ordine di grandezza di O(n3 /3).
13
Definizione 1.1.8 (Norma di un vettore). La norma ρ di un vettore v ∈ Rn è
definita come: (∑ )1/ρ
n
∥v∥ρ = |vi |ρ
i=1
Axi = λi xi ∀i = 1, 2, . . . , n
Ax = λx
Ax − λIx = 0
(A − λI)x = 0
p(λ) ∈ Pn : det(A − λI) = 0
Da notare che gli autovalori di una matrice diagonale (Def. 1.1.1), triangolare
inferiore (Def. 1.1.2) o triangolare superiore (Def. 1.1.3) sono gli elementi presenti
sulla diagonale principale.
Inoltre per una generica matrice A ∈ Rn×n invertibile vale la seguente proprietà:
1
λi (A−1 ) = ∀i = 1, 2, . . . , n (1.8)
λn+1−i (A)
14
Teorema 1.1.3 (Autovalori per matrice definita positiva). Data una matrice
A ∈ Rn×n simmetrica e definita positiva, i suoi autovalori sono tutti numeri reali
positivi:
{λi (A)}ni=1 ⊂ R+
Dalla Def. 1.1.11 si può ricavare che data la generica matrice A ∈ Rn×n :
√
∥A∥2 = λmax (AT A)
KP = ∥A∥P A−1 P
P ∈ [0, +∞)
15
Tabella 1.1: Principali stime per l’errore sulla soluzione di sistemi lineari.
Errore assoluto e = x − x̂ e ∈ Rn
∥e∥
Errore relativo er = ∥x∥
er ∈ R, x ̸= 0
Residuo r = b − Ax̂ r ∈ Rn
∥r∥
Residuo normalizzato rr = ∥b∥
rr ∈ R, b ̸= 0
Ora che sono state definite tutte i parametri utili, si può dire che:
Ax = b −→ (A + δA)x̂ = b + δb
16
Nonostante tutto è sempre possibile avere una stima di massima dell’errore, in
funzione del numero di condizionamento (Def. 1.1.12) dalla matrice A:
[ ]
1 ∥δA∥2 ∥δb∥2
er ≤ +
∥δA∥
1 − K2 (A) ∥A∥ 2 ∥A∥2 ∥b∥2
2
Quindi se non ci sono fluttuazioni dei coefficienti della matrice del sistema lineare
δA = 0 si ottiene:
er ≤ K2 (A)rr
Da questo nasce il principio per cui per sistemi lineari mal condizionati, caratte-
rizzati quindi da numeri di condizionamento K2 (A) alti, anche a fronte di residui
piccoli, si hanno grandi errori.
17
dove B è la matrice di iterazione e g il vettore di iterazione. Affinché il metodo
iterativo converga alla soluzione deve valere l’ipotesi di consistenza forte (Def. 3 a
pagina XII), che implica che devono valere le seguenti equazioni:
x = A−1 b x = Bx + g
g = x − Bx = (I − B)x = (I − B)A−1 b
Per avere la convergenza alla soluzione l’errore si deve ridurre con l’aumentare
delle iterazioni:
e(k) ≤ ∥Bk ∥2 ∥e(0) ∥ per k = 1, 2, . . . (1.10)
18
In più, tanto più e piccolo il raggio spettrale ρ(B) → 0, tanto più è rapida la
convergenza.
aii ̸= 0 ∀i = 1, 2, . . . , n
Di conseguenza si ha che:
x(0) ∈ Rn arbitrario
Per k = 1, 2, 3 . . .
Per i = 1, 2, . . . , n
[ ]
(k+1) 1 ∑n
(k)
xi = bi − aij xj
aii j=1
j̸=i
19
Teorema 1.2.1 (Convergenza sufficiente al metodo di Jacobi). La condizione
sufficiente (ma non necessaria) per risolvere un sistema con il metodo di Jacobi è
che la matrice A ∈ Rn×n sia a dominanza diagonale stretta per righe (Def. 1.1.5 a
pagina 7) oppure simmetrica e definita positiva.
A=D−E−F
PGS = D − E
BGS = I − (D − E)−1 A
20
Infine, l’algoritmo di Gauss-Seidel è:
x(0) ∈ Rn arbitrario
Per k = 1, 2, 3 . . .
Per i = 1, 2, . . . , n
[ ]
(k+1) 1 ∑
i−1
(k+1)
∑n
(k)
xi = bi − aij xj − aij xj
aii j=1 j=i+1
Importantissimo notare che l’algoritmo non è in alcun modo implicito, poiché tutte
(k+1)
le componenti xj sono già state calcolate.
Esattamente come per il metodo di Jacobi, si ha che la condizione sufficiente,
ma non necessaria è avere una matrice A o simmetrica e definita positiva, oppure
a dominanza diagonale stretta per righe.
• simmetrica A = AT ;
21
Teorema 1.2.2 (Relazione metodo del gradiente con l’energia del sistema). Il
vettore x è soluzione del sistema lineare Ax = b se e solo se è punto di minimo
assoluto dell’energia del sistema ϕ(y).
1
ϕ(y) = (x + e)T A(x + e) − (x + e)T b
2
1 T 1 1 1
= x Ax + xT Ae + eT Ax + eT Ae − xT b − eT b
2 2 2 2
1 T 1 T 1 T 1
= x Ax − x b + e Ae + x Ae + eT Ax − eT b
T
|2 {z } |2 {z } |2 {z } 2
=ϕ(x) = 12 ∥e∥A k∈R⇒k=kT
( )T
1 1 T 1
= ϕ(x) + ∥e∥A + x Ae + eT Ax − eT b
2 2 2
1 1 1
= ϕ(x) + ∥e∥A + eT |{z}
AT x + eT Ax − eT b
2 2 T
2
A =A
1 1 1
= ϕ(x) + ∥e∥A + eT Ax + eT Ax − eT b
2 2 2
1 ( )
= ϕ(x) + ∥e∥A + eT Ax − b
2 | {z }
=0 (soluz.)
1
= ϕ(x) + ∥e∥A ≥ ϕ(x)
2
dato che x è minimo quindi, per il Teorema di Fermat, deve anche essere un punto
22
stazionario; questo implica che la funzione ha gradiente nullo nel punto x:
∇ϕ(y) x = 0
A sua volta questo implica che tutte le componenti di ϕ(y) si annullano nel punto.
Riscrivendo poi la funzione descitta dall’Eq. (1.11) come somma di n componetenti
e derivando la generica k-esima nel punto si ottiene:
1∑ ∑
n n
1 T
ϕ(y) = y Ay − y b =
T
yi aij yj − y i b1
2 2 i,j=1 i=1
1∑ 1∑
n n
∂ϕ(y)
= akj yj + yi aik − bk
∂yk 2 j=1 2 i=1
∑
n
= akj yj − bk
j=1
23
a questo punto serve calcolare l’incremento αk ad ogni iterazione: per farlo si può
utilizzare una funzione scalare F (α) : R −→ R definita come:
( )
F (αk ) = ϕ(x(k+1) ) = ϕ x(k) − αk ∇ϕ(x(k) )
Derivando la funzione F (α) e tenendo presente che per definizione (Eq. (1.13))
∇ϕ(x(k) ) = −r(k) si ottiene il prodotto scalare di:
( )
′
( (k) ) T
F (αk ) = − ∇ϕ x − αk ∇ϕ(x ) (k)
∇ϕ(x(k) )
( )
( (k) )
= −∇ϕ(x ) ∇ϕ x − αk ∇ϕ(x )
(k) T (k)
T
= r(k) ∇ϕ(x(k) − αk r(k) )
Volendo infine avere αk come punto di ottimo, si impone che la derivata prima
F ′ (αk ) sia nulla (Teo. Fermat), da questo ricava l’incremento αk :
T T
F ′ (αk ) = −r(k) r(k) + r(k) Ar(k) = 0
T
r(k) r(k)
αk = T
r(k) Ar(k)
A questo punto si hanno tutti gli elementi necessari per scrivere l’algoritmo del
24
metodo del gradiente semplice:
x(0) arbitrario
r(0) = b − Ax(0)
Per k = 1, 2, . . .
T
r(k) r(k)
αk = T
r(k) Ar(k)
(k+1)
x = x(k) + αk r(k)
r(k+1) = b − Ax(k+1)
Teorema 1.2.3 (Errore metodo del gradiente). Sia A una matrice simmetrica e
definita positiva. Si ha che il metodo del gradiente converge a alla soluzione x
∀x ∈ Rn e che:
K(A) − 1
e(k) A
≤ dkp ∥e(0) ∥A con dp =
K(A) + 1
25
Dove p(k) serve ad individuare la direzione di massima decrescita in modo tale da
garantire la soluzione nel minor numero di passaggi possibili, mentre αk serve a
garantire che ϕ(x(x+1) ) sia punto di minimo. Come per il metodo del gradiente,
αk minimizza la funzione ϕ(x(x+1) ) solo quando:
T
(k) p(k) r(k)
α = T
p(k) Ap(k)
pk Apj = 0 ∀j = 1, 2, . . . , k − 1 ∀k = 2, 3, . . . , n
dove α1 e α0 ancora una volta tali che la funzione descritta dall’Eq. (1.11)
sia minima:
ϕ(x(0) + α0 p(0) + α1 p(1) ) = F (α0 , α1 )
26
dove F è una funzione da due a una sola variabile.
4. Come nel caso precedente, si deve andare trovare il minimo della funzione
ϕ(y), che, a differenza del caso precedente, dipende da due parametri distinti:
α0 e α1 . Per calcolarli serve per prima cosa andare ad azzerare il gradente
∇F (α0 , α1 ):
∂F T
= p(0) ∇ϕ(x(0) + α0 p(0) + α1 p(1) )
∂α0
T[ ( ) ]
= p(0) A x(0) + α0 p(0) + α1 p(1) − b
T ( )
= p(0) −r(0) + α0 Ap(0) + α1 Ap(1)
T T T
= −p(0) r(0) + α0 p(0) Ap(0) + α1 p(0) Ap(1) = 0
∂F T T T
= · · · = −p(1) r(0) + α0 p(1) Ap(0) + α1 p(1) Ap(1) = 0
∂α1
T
p(0) Ap(1) = 0
tale condizione implica che la direzione p(1) sia A-coniugata (Def. 1.2.2)
rispetto alla precedente. Infine, si può calcolare l’ultimo parametro incognito
α1 :
T
p(1) r(1)
α1 = T
p(1) Ap(1)
27
Nel caso generale il principio non cambia:
∑
k
{ }
(k+1) (0)
x =x + αj p(j) = min ϕ(y(k+1) ) = F (α)
y(k+1) ∈Rn
j=0
p(k)T Ap(j) = 0 ∀j = 1, 2, . . . , k − 1
αk tale che
∇F (α) = 0
T
p(k) r(k)
=⇒ αk = T
p(k) Ap(k)
A questo punto non rimane che trovare i p(k) A-coniugati tra loro: per fare questo,
alla fine di ogni passaggio si modifica la direzione di ∇ϕ (che corrisponde all’op-
posto della direzione del residuo) sottraendo un multiplo del vettore p(k) , ovvero
si cerca p(k+1) della forma:
Ovviamente si sceglie βk tale che p(k+1) sia A-coniugato a tutti i {p(j) }kj=1 prece-
denti. Di fatto non si deve risolvere un sistema di k incognite, ma si βk ricava da
un’unica equazione scalare:
T T ( )
p(k) Ap(k+1) = p(k) A r(k+1) − βk p(k) = 0
T
p(k) Ar(k+1)
=⇒ βk = T
p(k) Ap(k)
x(0) arbitrario
r(0) = b − Ax(0)
p(0) = r(0)
28
Per k = 1, 2, . . .
T
r(k) p(k)
αk = T
p(k) Ap(k)
x(k+1) = x(k) + αk p(k)
r(k+1) = b − Ax(k+1)
T
p(k) Ar(k)
βk = T
p(k) Ap(k)
p(k+1) = r(k+1) − βk p(k)
pendenti fra loro. Questo fatto a sua volta implica che ogni volta il metodo calcola
una nuova direzione e la utilizza per minimizzare la funzione ϕ(y) (cfr. Eq. 1.11
a pagina 21) aggiunge una dimensione alla direzione di decrescita: quindi prima
minimizza con una retta, poi con un piano, poi uno spazio, poi un iperpiano di
quattro dimensioni. . . fino ad arrivare alla soluzione esatta in massimo n passaggi.
Per questo motivo si può dire che il metodo del gradiente coniugato è un metodo
diretto, ma la pratica comune accetta senza particolari problemi un’approssima-
zione della soluzione: ad esempio per risolvere un sistema di mille equazioni lineari
non sono necessarie mille iterazioni, ma solo il numero minimo necessario per rag-
giungere la tolleranza prefissata, soprattutto perché a un certo punto gli errori di
calcolo diventano preponderanti.
Teorema 1.2.4 (Errore metodo del gradiente coniugato). Sia A una matrice sim-
metrica e definita positiva. Si ha che il metodo del gradiente coniugato converge
a alla soluzione x ∀x(0) ∈ Rn in al più n iterazioni e che l’errore e(k) è ortogonale
a tutte le direzioni di decrescita e ha norma:
√
2ck K(A) − 1
e (k)
A
≤ 2k
∥e ∥A
(0)
con c = √ per k < n
1+c K(A) + 1
Anche in questo caso, visto che si tratta di un algoritmo molto potente, è già
implementato nell’ambiete MATLAB® , tramite la funzione:
29
x = pcg(A, b, tol)
dove A è la matrice del sistema, b il vettore dei termini noti, x la soluzione (output
della funzione) e tol la tolleranza richiesta.
Questa funzione è anche più potente di quella che è stata appena descritta,
perché fa uso della tecnica di precodizionamento per migliorare il numero di con-
dizionamento (Def. 1.1.12) della matrice A e quindi per facilitare la risoluzione del
sistema, evitando inutili errori computazionali.
La regola generale per gli stimatori di errore ẽ(k) afferma che uno stimatore è
efficace solo quando:
ẽ(k) ≃ e(k+1)
e(k+1) ≤ K2 (A)ẽ(k)
30
da cui segue che:
e(k) ≤ e(k+1) + δ (k)
31
Capitolo 2
All’interno del Cap. 1 sono stati analizzati vari metodi per risolvere i sistemi lineari.
Ma come si trova la soluzione di un’equazione o sistema di equazioni che non
godono della proprietà di linearità? In questo capitolo si analizzeranno i metodi
per rispondere a questa domanda.
Consideriamo un problema molto comune nel campo dell’ingegneria, cioè quello
relativo all’esistenza (e alla ricerca) degli zeri di una funzione a variabile reale
generica:
f (x) : [a, b] ⊆ R −→ R
α ∈ R | f (α) = 0
In genere questo problema rimane circoscritto a un intervallo specifico [a, b], sia
per motivi numerici, sia per motivi tecnico-ingegneristici: non è molto significativo
andare a trovare uno zero i un intervallo che risulta fuori-range per un determinato
fenomeno fisico. Prima di procedere con la trattazione ricordiamo un teorema
fondamentale per la risoluzione di equazioni.
Teorema 2.0.1 (Degli zeri). Sia f (x) ∈ C(I) con I = [a, b] ⊆ R tali che
f (a)f (b) < 0. Allora esiste almeno un punto interno all’intervallo I in cui la
funzione f (x) si annulla:
∃ α ∈ (a, b) | f (α) = 0
32
Nelle sezioni seguenti vengono discussi i principali metodi per risolvere le equa-
zioni implicite. Nell’ambiente MATLAB® è già implementata una funzione in
grado di risolvere per via numerica un’equazione:
x = fsolve(f, x0)
s = solve(eqn, var)
che permette di risolvere equazioni simboliche per via analitica: in questo caso
eqn è l’equazione da azzerare, mentre var è la variabile rispetto a cui ricavare la
soluzione s.
33
Se f (a(k) )f (x(k) ) > 0
a(k+1) = x(k)
b(k+1) = b(k)
a(k) + b(k)
x(k+1) =
2
Come si nota dal testo, questo metodo non fa altro che prendere l’intervallo
iniziale [a, b], valutare la funzione nel punto medio e scegliere quale tra i due nuovi
intervalli rispetta condizioni del teorema degli zeri (Teo. 2.0.1); il metodo si arresta
quando si giunge al di sotto di una tolleranza prefissata, oppure se a un certo punto
la funzione viene valutata proprio in corrispondenza del suo zero. In ogni caso, ad
ogni iterazione l’intervallo si va a dimezzare:
|I|
I (k) = Per k = 1, 2, . . . (2.1)
2k
Questa equazione è in qualche modo indice dell’errore massimo che il metodo sta
commettendo al passaggio k; quest’ultimo infatti non può in nessun caso essere
maggiore della metà della misura dell’intervallo I (k) :
I (k) |I|
e(k) = α − x(k) ≤ ẽ(k) = = k+1
2 2
|I| k→+∞ +
ẽ(k) = −−−−→ 0
2k+1
1
f (x) = √ + ln(x + 1) α = 1,020 747
x+1
Inoltre, sempre dell’Eq. (2.1), si può anche calcolare il numero minimo di passaggi
34
100
10−4
10−8
10−12
10−16
0 10 20 30 40 50
Figura 2.1: Esempio dell’andamento degli errori del metodo di bisezione in funzione
del numero delle iterazioni su un grafico semilogaritmico: in particolare, in azzurro è
riportato l’errore effettivo, mentre in arancione l’errore stimato.
b−a b−a
ẽ(k) = < toll =⇒ k > log2
2k+1 toll
x(k+1) − α
lim ( (k) )p = µ
k→+∞ x −α
In generale quindi, si cerca di avere p maggiore possibile, per avere una con-
vergenza più veloce allo zero della funzione. Nel caso del metodo di bisezione il
metodo converge con andamento lineare (p = 1), infatti:
35
2.2 Metodo di Newton
2.2.1 Metodo di Newton semplice
Il metodo di Newton è uno strumento potente per ricercare uno zero α di una
data funzione f ∈ C 1 (I) con I = [a, b] ⊆ R. Il principio di base è quello di
approssimare la soluzione con l’intersezione della retta tangente al un punto della
funzione con l’asse delle ascisse. Per ricavare il metodo, si consideri l’espansione
con il polinomio di Taylor dell’intorno di x0 arrestata al primo ordine di una
generica funzione f ∈ C 1 (I):
( )
y(x) = f (x0 ) + f ′ (x0 )(x − x0 ) + o (x − x0 ) (2.2)
f (x0 )
x = x0 −
f ′ (x0 )
36
vicino a α la convergenza del metodo di Newton è garantita.
Se inoltre f ∈ C 2 (I) il metodo di Newton converge con ordine p = 2:
x(k+1) − α 1 f ′′ (α)
lim ( ) =
k→+∞ x(k) − α 2 2 f ′ (α)
In particolare si dice che α è uno zero semplice se m = 1, mentre che α è uno zero
multiplo se m > 1. Inoltre la radice α ha molteplicità m se e solo se f (x) può
essere riscritta come:
x(k+1) − α 1
lim = 1 −
k→+∞ x(k) − α m
37
Il metodo di Newton può essere riscritto come una funzione ϕ(x) tale che:
f (x)
ϕ(x) = x −
f ′ (x)
(x − α)m g(x)
ϕ(x) = x −
m(x − α)m−1 g(x) + (x − α)m g ′ (x)
(x − α)g(x)
=x−
mg(x) + (x − α)g ′ (x)
′
[ ′
]
mg 2
(x) + (x − α)g (x) 1 + mg(x) + (x − α)g (x)
ϕ′ (x) = 1 − ( )2 +
mg(x) + (x − α)g ′ (x)
mg ′ (x) + g ′ (x) + (x − α)g ′′ (x)
+ (x − α)g(x) ( )2
mg(x) + (x − α)g ′ (x)
mg 2 (α) 1
ϕ′ (α) = 1 − 2 2
=1− con m > 1
m g (α) m
Quindi il Teo. 2.2.2 afferma che tanto più è grande m tanto più è lenta la con-
vergenza del metodo. Tale complicazione spesso non è accettabile, perché impiega
grande potenza computazionale per risolvere un’equazione.
Per fortuna, è possibile migliorare la convergenza introducendo una costante m
davanti al rapporto della funzione con la sua derivata. L’algoritmo diventa quindi:
38
metodo di Newton modificato converge con p = 2.
F(x) : Rn −→ Rn
α ∈ Rn | F(α) = 0
A questo punto si può pensare di andare a generalizzare i metodi visti nei para-
grafi precedenti per le funzioni a una sola variabile reale. Purtroppo però il metodo
più semplice, cioè quello di bisezione (Sez. 2.1) non può essere implementato per
le funzioni con dominio di dimensione superiore a quella di R, perché non si riesce
a dividerlo ogni volta in due parti uguali. Rimane quindi il metodo di Newton
(Sez. 2.2.1) la cui versione base può essere facilmente estesa per domini di Rn .
Definizione 2.3.1 (Matrice Jacobiana). Data una funzione F ∈ C 1 (Rn ), si defi-
nisce matrice Jacobiana la matrice le cui componenti sono pari a:
[ ] ∂fi ∀i = 1, 2, . . . , n
JF = (x)
ij ∂xj ∀j = 1, 2, . . . , n
39
Ancora una volta trascurando i termini di ordine superiore al primo, imponendo
che la funzione sia nulla in x e ipotizzando che det(JF ) ̸= 0 si ottiene un sistema
lineare con incognita δx:
JF (x0 ) (x − x0 ) = −F(x0 )
| {z }
=δx
Come già analizzato prima, per ottenere l’algoritmo del metodo basta reiterare il
processo:
x(0) arbitrario
Per k = 0, 1, 2, . . .
JF (x0 (k) )δx(k) = −F(x(k) ) si calcola δx(k)
x(k+1) = x(k) + δx(k)
Per determinare quale sia il metodo migliore per risolvere il sistema lineare, si
rimanda al Cap. 1 a pagina 2.
40
riportato in Fig. 2.2, il punto fisso di una funzione può essere interpretato come
intersezione con la bisettrice del primo e terzo quadrante.
La grande utilità delle iterazioni di punto fisso è rappresentata dal fatto che
possono essere utilizzate per la ricerca di zeri di alcune funzioni. Infatti:
ϕ(α) = α =⇒ ϕ(α) − α = 0
Quindi il punto fisso α della funzione coseno corrisponde allo zero della funzione:
f (x) = ϕ(x) − x
Teorema 2.4.1 (Convergenza globale iterazioni di punto fisso). Sia ϕ(x) una
funzione a variabili reali ϕ(x) : R −→ R.
1. Se ϕ ∈ C 0 ([a, b]) e ϕ(x) ∈ [a, b] ∀x ∈ [a, b] allora esiste almeno un punto fisso
della funzione ϕ in [a, b].
3. Se inoltre, ϕ ∈ C 1 ([a, b]) con |ϕ′ (x)| < 1 ∀x(0) ∈ [a, b] allora esiste un solo
punto fisso e l’algoritmo converge almeno linearmente. Inoltre:
x(k+1) − α
lim = ϕ′ (α)
k→+∞ x(k) − α
Teorema 2.4.2 (Valor medio di Lagrange). Sia f ∈ C 1 ([a, b]) allora esiste almeno
un ξ ∈ [a, b] tale che:
f (a) − f (b) = f ′ (ξ)(a − b)
41
1.5 y
1
α
0.5
x
0.5 1 1.5 2
Figura 2.2: Punto fisso della funzione y(x) = cos(x), ottenuto dall’intersezione della
funzione cos(x) con la bisettrice del primo e quarto quadrante (tratteggiata).
Teorema 2.4.3 (Convergenza locale di Ostrowski). Se ϕ ∈ C 1 ([a, b]) e |ϕ′ (x)| < 1,
assegnato un x(0) arbitrario abbastanza vicino ad α allora il metodo iterativo di
punto fisso converge con ordine p = 1:
x(k+1) − α
lim = ϕ′ (α)
k→+∞ x(k) − α
che per il teorema del valore medio di Lagrange (Teo. 2.4.2) è pari a:
x(k+1) − α
= ϕ′ (ξ (k) )
x(k) − α
Inoltre sapendo che per k → +∞, x(k) → α e ξ (k) → α, poiché il valore di ξ (k) è
42
limitato all’intervallo [x(k) , x(k+1) ]; quindi per il teorema del confronto semplice:
Quindi se |ϕ′ (x)| < 1 (una delle ipotesi del teorema) si avrà la convergenza, se
invece |ϕ′ (x)| > 1 il metodo non arriverà alla soluzione, a meno che x(0) ∼
= α. Vi
′
è poi un caso limite dato da |ϕ (x)| = 1, in questo caso la convergenza del metodo
dipende solamente dalla scelta della prima iterata x(0) .
x(k+1) − α 1 ′′
lim ( )2 = ϕ (α)
k→+∞ x (k) −α 2
α ∈ Rn | ϕ(α) = α
Anche i punti fissi di funzioni vettoriali possono essere utilizzati per ricercare lo
zero della funzione F(x) = ϕ(x) − x, tramite l’algoritmo:
x(0) arbitrario
Per k = 0, 1, 2, . . .
x(k+1) = ϕ(x(k) )
43
2.4.2 Metodo di Newton come iterazione di punto fisso
il metodo di Newton, già descritto nella Sez. 2.2.1 a pagina 36, può essere analizzato
come un caso particolare dei iterazione di punto fisso, infatti:
f (x(k) )
x(k+1) = x(k) − = ϕN (x(k) ) (2.5)
f ′ (x(k) )
f (x)f ′′ (x)
ϕ′N (x) = ( )2
f ′ (x)
x(k+1) − α 1 f ′′ (α)
lim ( )2 = con f ′ (α) ̸= 0
k→+∞ x(k) − α 2 f ′ (α)
f (x(k) )
x(k+1) = x(k) − m = ϕM (x(k) ) (2.6)
f ′ (x(k) )
44
2.5 Criteri di arresto
Come si è analizzato nei precedenti paragrafi, la maggior parte dei metodi numerici
per il calcolo degli zeri è di tipo iterativo, per cui si devono trovare dei parametri
adatti ad arrestare il processo quando il risultato ottenuto giunge al di sotto delle
tolleranze prefissate.
Per ovviare a questo problema esistono principalmente due diverse tipologie di
criteri:
Residuo a differenza di quanto visto per i sistemi lineari, questo criterio di arresto
è indice dell’errore assoluto che si sta commettendo (la soluzione è sempre
nota, poiché si ricerca sempre lo zero della funzione).
I criteri appena visti possono essere facilmente estesi alla risoluzione di sistemi
non lineari (Sez. 2.3 a pagina 39), sostituendo il valore assoluto, con la norma dello
zero.
45
Capitolo 3
∑
n
f (k) (x0 )
Tn (f, x) = (x − x0 )k
k=0
k!
46
e si ha che:
f (x) = Tn (f, x) + Rn (x)
f (n) (ξ)
Rn (x) = (x − x0 )n con ξ ∈ (x, x0 )
n!
taylortool
Oltre a questo, è anche disponibile una funzione per il calcolo del polinomio di
Taylor per funzioni a variabili simboliche:
47
e inoltre se è nota la funzione che definisce gli yi si può scrivere:
{ }n
(xi , f (xi )) i=0 xi ̸= xj ∀i ̸= j
f˜(xi ) = yi ∀i = 0, 1, . . . , n
• Polinomiali
f˜(x) = a0 + a1 x + a2 x2 + · · · + an xn
• Razionali
a0 + a1 x + a2 x 2 + · · · + an x n
f˜(x) =
b0 + b1 x + b2 x 2 + · · · + b m x m
• Trigonometriche
∑
n
f˜(x) = cj eijx
j=−n
• Esponenziali
f˜(x) = k1 ax1 + . . .
• Splines
Per semplicità si tratteranno solamente gli interpolanti che dipendono solo linear-
mente dai loro coefficienti incogniti; quindi dall’elenco precedente sono escluse le
sole interpolazioni razionali ed esponenziali.
48
particolare per i calcolo degli zeri, integrali, derivate. . . ), ma purtroppo si vedrà
che non sempre riescono a seguire bene la funzione che interpolano. Prima di
procedere con i metodi di interpolazione polinomiale si analizzano le proprietà
dell’insieme P.
pn (x) = a0 + a1 x + a2 x2 + · · · + an xn
pn (x) ∈ Pn
∃! Πn (x) ∈ P
Πn (xi ) = yi ∀i = 0, 1, . . . , n
Data inoltre una funzione f (x) ∈ C 0 (I) tale che f (xi ) = yi , esiste un solo
polinomio di grado al più n che interpola f (x) nell’intervallo I:
∃! Πn f (x) ∈ P Πn f (xi ) = yi ∀i = 0, 1, . . . , n
49
Dimostrazione 3.1.1 (Unicità polinomio interpolante). Si supponga per assurdo
che vi siano due diversi polinomi dello stesso grado (al più n) che interpolano gli
stessi punti:
∃ Π1n (x), Π2n (x) ∈ P con Π1n (x) ̸= Π2n (x)
Dato che entrambi i polinomi hanno grado al più n una loro generica combinazione
lineare sarà ancora un un polinomio di grado al più n; tra le infinite possibili si
scelga:
.
Πn (x) = Π1n (x) − Π2n (x) Πn (x) ∈ Pn
Ora per definizione di polinomio interpolante (presente all’interno del Teo. 3.1.1)
si ha che:
Naturalmente la base non è unica, ma, come per tutti gli spazi vettoriali, ne
esistono infinte diverse in grado di generare l’intero Pn .
50
definisce base di Lagrange o polinomio caratteristico di Lagrange:
∑
n
Πn (x) = yi φi (x)
i=0
1. Deve essere nullo in tutti i k nodi diversi dall’i-esimo nodo che deve inter-
polare, per cui:
∏
n
φi (x) = H (x − xk )
k=0
k̸=i
51
2. Deve valere uno nel’i-esimo nodo che deve interpolare; imponendo questa
condizione si ottiene:
∏
n
1
φi (xi ) = 1 =⇒ H (xi − xk ) = 1 =⇒ H = ∏n
k=0 (xi − xk )
k=0 k̸=i
k̸=i
Per queste ragioni il generico polinomio che interpola gli n + 1 nodi assegnati è
esprimibile come:
∏n
∑
n ∑
n k=0 (x − xk )
yi ∏n
k̸=i
Πn (x) = yi φi (x) =
i=0 i=0
k=0 (xi − xk )
k̸=i
(x − . . . )(x − . . . )
φa (x) = = ...
(a − . . . )(a − . . . )
(x − (−1))(x − 0) x2 + x
φ0 (x) = =
(1 − (−1))(1 − 0) 2
(x − 1)(x − 0) x −x
2
φ1 (x) = =
(−1 − 1)(−1 − 0) 2
(x − 1)(x − (−1))
φ3 (x) = = −x2 + 1
(0 − 1)(0 − (−1))
Una volta trovate la basi si moltiplica per l’ordinata del rispettivo punto si som-
52
3 y
2
P1
1
P2 P0
x
−2 −1 1 2
−1
mano tutte (notare che si poteva evitare di calcolare φ2 (x), poiché y2 = 0).
x2 + x x2 − x 3 1
Π2 (x) = 1 +2 + 0(−x2 + 1) = x2 − x
2 2 2 2
53
Teorema 3.1.2 (Errore di interpolazione). Dati {xi }ni=0 ⊂ I noti distinti, che
appartengono a una funzione f ∈ C n+1 (I) allora l’errore di interpolazione è pari
a: ( )
f (n+1) ξ(x)
En f (x) = ωn (x)
(n + 1)!
con ξ(x) ∈ I. Il termine ωx (x) ci chiama polinomio nodale ed è pari a:
∏
n
ωn (x) = (x − xi )
i=0
1
en (f ) ≤ max f (n+1) (x) max|ωn (x)|
(n + 1)! x∈I x∈I
La formula che propone il Teo. 3.1.2 è molto generica e in generale poco utile,
serve quindi riarrangiarla in qualche modo per ottenere una grandezza utile a
risolvere i problemi di calcolo. Per farlo però, bisogna definire la tipologia dei nodi
che si vuole utilizzare.
Nodi equispaziati
Allora il valore massimo (in modulo) che può assumere il polinomio nodale ωn (x)
è pari a:
n!
max|ωn (x)| = hn+1
x∈I 4
Quindi il massimo valore dell’errore di interpolazione vale:
hn+1
en (f ) = max f (n+1) (ξ)
4(n + 1) ξ∈I
54
Dimostrazione 3.1.3 (Errore interpolazione con nodi equispaziati). Il problema
principale è sistemare il valore del polinomio nodale, in modo da poter trovare una
limitazione superiore al valore. Si consideri dunque il suo valore assoluto:
∏
n ∏
n
| ωn (x)| = (x − xi ) = |x − xi |
i=0 i=0
Supponga che il generico valore x appartenga al primo intervallo dei nodi [x0 , x1 ].
Si sa per certo che se x = x0 , x1 la moltiplicazione dei primi due termini della
produttoria sarà nulla. Inoltre, essendo una moltiplicazione del valore assoluto di
due polinomi la funzione sarà continua. Quindi la funzione |x − x0 ||x − x1 | avrà
(per il teorema di Rolle) almeno un punto di massimo, che per ovvie ragioni sarà
al centro dell’intervallo. In formule:
hh h2
|x − x0 ||x − x1 | ≤ =
22 4
Dato che x ∈ [x0 , x1 ] le distante da tutti gli altri punti sono facilmente calcolabili:
|x − x2 | ≤ 2h |x − x3 | ≤ 3h ... |x − xn | ≤ nh
∏n
h2 hn+1
| ωn (x)| = |x − xi | ≤ n! hn−1 = n!
i=0
4 4
Il Teo. 3.1.3 e più in generale il Teo. 3.1.2 non garantiscono in alcun modo che
aumentando il numero dei nodi l’errore dell’interpolazione diminuisca, perché il
massimo della derivata (n + 1)-esima potrebbe tendere a infinito più velocemente
55
2 y
0.6 y
0.4
1
0.2
x
x
−3 −2 −1 1 2 3
−3 −2 −1 1 2 3
−0.2
−1
(a) Funzione di Runge (azzurro) con con po- (b) Andamento dell’errore En f (x): in aran-
linomi di Lagrange di grado 5 (arancione) e 7 cione per il polinomio di Larange di grado 5
(magenta). e in magenta per quello di grado 7
[ ]−1
di quanto ωn (x) (n + 1)! tenda a zero. La regola dice che se:
n→+∞
max f (n+1 è limitato =⇒ en (f ) −−−−→ 0+
x∈I
Fenomeno di Runge
Esistono delle funzioni, anche apparentemente semplici, tali per però sono sempre
affette da un grande errore di interpolazione. Un esempio di esse è la funzione di
Runge, definita come:
1
f (x) =
1 + x2
Nonostante non abbia una forma particolarmente complessa da studiare, non si
presta per nulla bene all’interpolazione di Lagrange, sopratutto sugli estremi. Un
esempio di questo fenomeno è riportato in Fig. 3.2.
56
Si osserva che anche per piccolissime variazioni l’interpolazione si allontana di
molto dalla funzione originale, nel caso più generico possibile si scrive:
Questi possono essere estesi al generico intervallo [a, b] dilatandoli per la metà della
misura dell’intervallo e translandoli di un valore pari alla media aritmetica degli
estremi: ( )
a+b b−a πi
x1 = − cos ∀i = 0, 1 . . . , n
2 2 n
Teorema 3.1.5 (Errore interpolazione per nodi CGL.). Se la funzione da in-
terpolare f ∈ C (n+1) (I) e si usano i nodi CGL. il fenomeno di Runge non si
verifica:
lim en (f (x)) = 0
n→+∞
57
Tornando al problema della costruzione del polinomio interpolante si può pen-
sare di ricavare Πn (x) vincolandolo a passare per tutti gli n + 1 punti assegnati;
questo significa risolvere un sistema lineare del tipo:
a0 + a1 x + a2 x2 + · · · + an xn = y0
a + a x + a x 2 + · · · + a x n = y
0 1 2 n 1
.
.. .
.. .
..
a + a x + a x 2 + · · · + a x n = y
0 1 2 n n
=⇒ Ba = y
dove a è il vettore dei coefficienti incogniti, y le ordinate dei punti nodali e infine
B la matrice di Vandermonde, cioè una matrice della forma:
1 x1 x21 x31 . . . xn1
1 x2 x22 x32 . . . xn2
V=
.. .. ..
. . .
1 xn xn x3n . . .
2
xnn
Ancora si può avere una prova che il polinomio di Lagrange esiste ed è unico,
poiché il determinante della matrice B può essere scritto come:
∏
det(B) = (xi − xj )
1≤i<j≤n
ed essendo tutti i nodi distinti per ipotesi det(B) ̸= 0 e quindi per il teorema
di Cramer esiste una sola soluzione al sistema e di conseguenza al problema di
interpolazione. Il problema di fondo è che la matrice di Vandermonde è mal con-
dizionata (K2 (B) ≫ 0) e, come visto nel Cap. 1 questo implica che i coefficienti
trovati saranno affetti da notevoli errori di calcolo.
58
costituita dal fatto che in questo caso l’intervallo I = [a, b] viene diviso in tanti
[Link] del tipo:
Ii = [xi−1 − x1 ] ∀ii = 1, 2, . . . , n
È importante notare che non è stata fatta nessuna ipotesi in merito all’ampiezza
Hi degli intervalli. Di conseguenza si definisce ampiezza caratteristica l’ampiezza
dell’intervallino di misura massima:
H = max |Hi |
i=1,2,...,n
1 f (x) − f (x)
e1H (f ) = max ΠH
x∈I
H2
1 (f ) ≤ ẽ1 (f ) =
eH H
max|f ′′ (x)|
8 x∈I
59
y
10
x
−2 2 4 6
−10
Figura 3.3: Interpolazione lineare composita della funzione f (x) = 3x sin(x) − e−x ,
riportata in azzurro: nello specifico la spezzata arancione è stata costruita utilizzando
cinque nodi equispaziati, quella magenta utilizzandone nove.
r (x) ∈ Pr
ΠH con x ∈ Ii ∀i = 1, 2, . . . , n
60
zione f ∈ C r+1 (I) allora l’errore di interpolazione è pari a:
1 (f ) ≤ ẽ1 (f ) = Cr H
eH con Cr ∈ R+
H r+1
max f (r+1) (x)
x∈I
Sr (x) Ii
∈ Pr ∀i = 1, 2, . . . , n
Sr (xj ) = yj ∀i = 1, 2, . . . , n
(k) − ∀i = 1, 2, . . . , n
Sr(k) (x+
j ) = Sr (xj )
∀k = 1, 2, . . . , k
In generale esistono molte tipologie diverse di spline, come quella naturale, che
prevede che la derivata seconda del primo e dell’ultimo nodo sia nulla o quella
not-a-knot, la quale richiede che tutte le derivate terze di ogni nodo siano uguali
a destra e a sinistra. MATLAB® a meno di diverse istruzioni utilizza l’ultima ed
è richiamabile attraverso il comando:
PP = spline(x, y)
61
In cui il vettore PP contiene i coefficienti del polinomio, mentre x e y sono le
coordinate dei punti nodali. Se si vuole avere una function handle che dato un
qualsiasi valore appartenente all’intervallo restituisca il valore dell’interpolazione
di tipo spline basta inserire su MATLAB® l’espressione:
s = @(t) spline(x, y, t)
A questo punto basta richiamarla nel codice tramite s(t), per ottenere il valore
ricercato.
∑
n
[ ]2 ∑
n
[ ]2
˜
yi − fm (xi ) ≤ yi − pm (xi ) ∀pm (x) ∈ Pm
i=0 i=0
62
0
−10
−20
−30
0 20 40 60 80 100
Figura 3.4: Il grafico mostra un esempio di un’interpolazione lineare nel senso dei
minimi quadrati di cento punti.
Definizione 3.2.2 (Energia del vettore dei coefficienti). Dato il vettore dei coef-
ficienti b ∈ Rm si chiama energia del vettore la funzione:
ϕ(b) : Rm −→ R
∑
n
[ ]2
ϕ(b) = yi − (b0 + b1 x1 + · · · + bm xm
i ) (3.1)
i=0
Teorema 3.2.1 (Calcolo polinomio approssimante nel senso dei minimi quadrati).
Il polinomio approssimante nel senso dei minimi quadrati di un insieme di n+1 dati,
con x1 sono nodi distinti, è un polinomio di grado m il cui vettore dei coefficienti
minimizza la funzione ϕ(b) definita dall’Eq.(3.1):
[ ]
a ∈ Rm | ϕ(a) = minm ϕ(b)
b∈R
63
te un punto stazionario, questo sarà un minimo assoluto. Per trovare i coefficienti
a si deve risolvere l’equazione vettoriale:
∇ϕ(b) a
=0 (3.2)
Il polinomio può essere riscritto come somma di termini che vanno da zero (termine
noto) a m:
∑n [ ∑ m ]2
∂ϕ ∂
(b) = yi − j
bj x i
∂bk i=0
∂bk j=0
∑ n [ ∑m ]
=2 yi − j
bj xi xki
i=0 j=0
∑
n ∑
m ∑
n
aj xij+k = xki yi ∀k = 0, 1 . . . , m
i=0 j=0 i=0
Aa = q
a = (a0 , a1 , . . . , am )T a ∈ Rm (incognito)
∑n
[A]kl = xk+l
i ∀k, l = 0, 1, . . . , m A ∈ R(m+1)×(m+1)
i=0
∑
n
qk = xki yi ∀k = 0, 1, . . . , m q ∈ Rm
i=0
64
Importante notare che tutte sommatorie sull’indice i sono indipendenti dalla di-
mensione del problema lineare. Nello specifico l’indice i va da zero a n ed è relativo
ai nodi, gli indici k e l per la matrice e k per il vettore dei termini noti sono relativi
al sistema lineare. Per fare maggiore chiarezza si può riscrivere il sistema come:
∑ ∑ ∑
n+1 xi ...
m
i xi
a y
∑ ∑i 2
1
∑ i
i
i xi
a1
. i xi i y x
i i
. .. .. .. = ..
. . . . .
∑ m ∑ 2m ∑ m
i xi ... x
i i a m i y i x i
Si nota molto bene che la matrice A è simmetrica e definita positiva (se tutti i
nodi xi sono distinti) quindi il sistema è sempre risolvibile.
In merito agli errori commessi dal metodo dei minimi quadrati non si fanno
grandi considerazioni. In genere comunque si calcola il intervallo interpolante
tale per cui un possibile risultato di un esperimento ricada all’interno di esso con
probabilità nota.
Nell’ambiente MATLAB® sono invece disponibili due funzioni utili per inter-
polare un insieme di punti. La prima consente di trovare i coefficienti coeff un
polinomio di grado n che interpoli nel senso dei minimi quadrati un insieme di
punti di coordinate x, y:
coeff = polyfit(x, y, n)
z = polyval(p, x)
Infine se si vuole creare una function handle che, assegnato un valore t resituisca
il valore interpolato di un set di punti:
65
z = @(t) polyval(polyfit(x, f(x), n), t)
Sul sito di MATLAB® sono riportate con precisione tutte le altre funzionalità di
queste due applicazioni, come ad esempio il calcolo della media o della deviazione
standard dell’approssimazione ai minimi quadrati. . .
66
Capitolo 4
Integrazione e derivazione
numerica
67
ed, inoltre, è possibile, attraverso le proprietà della funzione integrale, creare un
algoritmo che restituisca il valore della primitiva di una funzione da un punto di
riferimento prefissato a un generico valore x.
∫ b ∫ x
I(f ) = f (t)dt =⇒ F (x) = f (t)dt
a x0
Iq (f ) ≃ I S (f˜) ∀x ∈ [a, b]
68
H su cui è definita la formula.
IPS M (f ) ≃ I(Π0 f )
Da cui segue che la formula di quadratura del punto medio semplice è definita
dalla seguente espressione:
a+b
IPS M (f ) = (b − a)f (x̄) con x̄ =
2
b−a
xk = a + kH ∀k = 0, 1, . . . , n H=
n
∑
n
xk + xk−1
IPCM (f ) = I(ΠH
0 f) = H f (x̄) con x̄k =
k=i
2
Teorema 4.1.1 (Errore formula del punto medio semplice). Se la funzione inte-
granda f ∈ C 2 ([a, b]) si ha che l’errore della formula del punto medio semplice è
pari a:
(b − a)3 ′′
eSP M (f ) = I(f ) − IPS M (f ) = f (ξ) ξ ∈ [a, b]
24
69
Teorema 4.1.2 (Valor medio integrale). Sia f : [a, b] ⊆ R −→ R integrabile.
Allora esiste un punto c ∈ [a, b] tale che:
∫ b
1
f (c) = f (x)dx
b−a a
Teorema 4.1.3 (Valor medio discreto). Questo teorema è una sorta di discretiz-
zazione del teorema precedente. Data una funzione f : [a, b] ⊆ R −→ R, esiste
almeno un punto c ∈ [a, b] tale che:
1 ∑
N
f (c) = f (xi )
b − a i=0
Dimostrazione 4.1.1 (Errore formula del punto medio semplice). L’errore della
formula del punto medio semplice può essere visto come la differenza tra l’in-
tegrale della funzione e l’integrale della funzione approssimata, in questo caso un
polinomio di grado zero, che può essere riscritta sfruttando le proprietà di linearità
dell’integrale:
∫ ∫
b b [ ]
eSP M (f ) = I(f ) − IPS M (f ) = f (x)dx − (b − a)f (x̄) = f (x) − f (x̄) dx
a a
′ f ′′ (ξ)
f (x) = f (x̄) + f (x̄)(x − x̄) + (x − x̄)2 ξ ∈ [a, b]
2
70
Sostituendo l’espressione ottenuta all’interno dell’integrale si ottiene:
∫ b [ ]
eSP M (f ) = f (x) − f (x̄) dx
a
∫
[ b
f ′′ (ξ) ]
= f (x̄) + f ′ (x̄)(x − x̄) + (x − x̄)2 − f (x̄) dx
a 2
∫ b
[ ′ f ′′ (ξ) ]
= f (x̄)(x − x̄) + (x − x̄)2 dx
a 2
∫ b ∫
′ f ′′ (ξ) b
= f (x̄) (x − x̄)dx + (x − x̄)2 dx
2
| a {z } a
=0 per il Teo. 4.1.2
[( )3 ( )3 ]
f ′′ (ξ) a+b a+b
= b− − a−
6 2 2
f ′′ (ξ) [ ]
= (2b − a − b)3 − (2a − a − b)3
48
(b − a)3 ′′
= f (ξ) ξ ∈ [a, b]
24
Teorema 4.1.4 (Errore formula del punto medio composita). Se la funzione in-
tegranda f ∈ C 2 ([a, b]) si ha che l’errore della formula del punto medio composita
è pari a:
(b − a)3 ′′
eC
PM (f ) = I(f ) − I C
PM (f ) = f (ζ) ζ ∈ [a, b]
24
Dimostrazione 4.1.2 (Errore formula del punto medio composita). L’errore della
formula del punto medio semplice può, come prima, essere riscritto sfruttando le
proprietà di linearità dell’integrale e della sommatoria:
∫ ∑ n ∫
∑
b n xk [ ]
eC
P M (f ) = I(f ) − IPCM (f ) = f (x)dx − H f (x̄k ) = f (x) − f (x̄k ) dx
a k=1 k=1 xk−1
∑
n
H3
eC
P M (f ) = f ′′ (ξk ) ξk ∈ [xk−1 , xk ]
k=1
24
71
Infine la sommatoria può essere riscritta come:
b − a 2 ′′
eC
P M (f ) = H f (ζ) ζ ∈ [a, b]
24
Infine, l’ordine di convergenza della formula del punto medio composita è pari a
due, infatti l’errore va a zero come il termine H 2
Come si può facilmente notare dal Teo. 4.1.1 e dal Teo. 4.1.4 gli errori del-
la formula di quadratura del punto medio semplice e composita sono nulli se la
derivata seconda della funzione f (x) è nulla in tutto l’intervallo di integrazione:
questo significa che la formula del punto medio ha grado di esattezza uno, perché
riesce a integrare in modo esatto tutte le costanti e le rette.
ITS (f ) ≃ I(Π1 f )
b − a( )
ITS (f ) = f (b) + f (a)
2
È immediato notare che la base del trapezio è data dalla differenza b − a, mentre
le due altezze dalla valutazione delle funzione negli estremi.
Nel caso composito l’area viene approssimata da n trapezi o, equivalentemente,
la funzione viene interpolata da un polinomio lineare composito:
ITC (f ) ≃ I(ΠH
1 f)
72
Nel caso di n + 1 nodi equispaziati la formula diventa:
∑
n
f (xk ) + f (xk−1 )
ITC (f ) =H
k=1
2
1 ( ) ∑
n−1
b−a
= H f (x0 ) + f (xn ) + H f (xk ) con H =
2 k=1
n
Teorema 4.1.5 (Errori formula di quadratura del trapezio). Data una funzione
f ∈ C 2 ([a, b]) allora l’errore della formula di quadratura del trapezio semplice vale:
(b − a)3 ′′
eST (f ) = I(f ) − ITS (f ) = − f (ξ) ξ ∈ [a, b]
12
In cui il segno negativo è un refuso dei calcoli di poca importanza (infatti sparisce
se si analizza il valore assoluto dell’errore). Per quella composita:
b − a 2 ′′
T (f ) = I(f ) − IT (f ) = −
eC ζ ∈ [a, b]
C
H f (ζ)
12
Quindi anche nel caso dei trapezi l’errore converge con ordine due.
Il Teo. 4.1.5 afferma che pur avendo utilizzato più informazioni e avendo speso
più potenza computazionale per risolvere l’integrale con la formula dei trapezi il
grado di esattezza è sempre uno e l’errore non solo non ha ordine di convergenza
migliore rispetto a prima, ma addirittura è maggiore di quello del punto medio!
73
e il punto medio dell’intervallo [a, b]. La formula che si ottiene è quindi:
b − a[ ] b+a
ISS (f ) = f (a) + f (b) + 4f (x̄) x̄ =
6 2
Nel caso della formula composita vengono scelti gli estremi di ogni sotto-intervallo
e il punto medio di ciascuno di essi. Nel caso di nodi equidistanti si ottiene:
H ∑[ ]
n
ISC (f ) = f (xk−1 ) + f (xk ) + 4f (x̄k )
6 k=1
xk + xk−1 b−a
x̄k = H=
2 n
(b − a)5 (4)
S (f ) = I(f ) − IS (f ) = −
eC ξ ∈ [a, b]
C
f (ξ)
2880
b − a 4 (4)
C (f ) = I(f ) − IC (f ) = −
eC ζ ∈ [a, b]
C
H f (ζ)
2880
Questa volta il grado di esattezza della formula è pari a tre, mentre l’ordine di
convergenza dell’errore è quattro. Inoltre la costante a denominatore fa si che gli
errori della formula di Cavalieri-Simpson siano notevolmente più piccoli rispetto
alle altre due formule di quadratura.
74
Tabella 4.1: Esempio di tabaulazione dei pesi e nodi di quadratura dalla formula di
Cavalieri-Simpson.
y0 = a α0 = (b − a)/6
y1 = (a + b)/2 α1 = 2(b − a)/3
y2 = b α2 = (b − a)/6
∑
n
Iq (f ) = I(f˜) = αj f (yj )
j=0
dove l’insieme degli {αj }nj=0 raccoglie tutti gli n + 1 pesi di quadratura, che sono
figli della funzione interpolante utilizzata, mentre {yj }nj=0 sono i nodi di quadratura
necessari alla formula.
Grazie a questa tipologia di formulazione è possibile tabulare i valori ᾱj e ȳj
relativi metodo utilizzato e calcolati sull’intervallo standard [−1, 1]. Durante l’ese-
cuzione dell’algoritmo questi vengono richiamati e adattati all’intervallo generico
[a, b].
b−a b−a a+b
αj = ᾱj yj = ȳj +
2 2 2
Nella Tab. 4.1 sono riportati, a titolo di esempio, pesi e nodi della formula di
Cavalieri-Simpson.
75
Tabella 4.2: Polinomi di Legendre per n ≤ 6.
Grado Polinomio
0 L0 (x) = 1
1 L1 (x) = x
2 L2 (x) = 12 (3x2 − 1)
3 L3 (x) = 12 (5x3 − 3x)
4 L4 (x) = 18 (35x4 − 30x2 − 3)
5 L5 (x) = 18 (63x5 − 70x3 − 15x)
6 L6 (x) = 1
16
(231x6 − 315x4 + 105x2 − 5)
Formule di Gauss-Legendre
L0 (x) = 1 L1 (x) = x
2k + 1 k
Lk+1 (x) = Lk (x) − Lk−1 (x)
k+1 k+1
I primi sette polinomi di Legendre sono riportati nella Tab.4.2 e sono graficati
nella Fig. 4.1.
76
y
1
x
−1.5 −1 −0.5 0.5 1 1.5
−1
Si può dimostrare che questi polinomi sono tutti ortogonali sull’intervallo [−1, 1]
e che sono una base dello spazio vettoriale Pn . Inoltre questi polinomi sono di gran-
de utilità nel calcolo numerico degli integrali. Infatti si può scrivere la formula di
quadratura gaussiana associata come:
∑
n
InGL (f ) = ᾱjGL f (ȳjGL )
j=0
in cui le quantità {ȳjGL }nj=0 , da adattare all’intervallo [a, b], sono gli zeri del
polinomio di Legendre di grado n + 1, ossia:
.
ȳjGL = x̄ | Ln+1 (x̄) = 0 ∀j = 0, 1, . . . , n
2
ᾱjGL = [ ][ ′ ]2 ∀j = 0, 1, . . . , n
1− (ȳjGL )2 Ln+1 (ȳjGL )
77
Tabella 4.3: Nodi, pesi e grado di esattezza si alcune formule di quadratura di Gauss-
Legendre. I pesi corrispondenti a nodi simmetrici rispetto allo zero sono riportati una
sola volta.
n Nodi Pesi r
0 0 2 1
1 ± √13 1 3
√
2 0, ± 515 5 8
,
9 9
5
.. .. .. ..
. . . .
Formule di Gauss-Legendre-Lobatto
Queste formule sono analoghe alle formule di Gauss-Legendre (Sez. 4.1.5), con la
sola differenza che il primo e l’ultimo nodo hanno valore fissato, indipendentemente
dal valore di n:
ȳ0 = −1 ȳn = −1
Questo vincolo ovviamente porta benefici di calcolo, al netto di una riduzione del
grado di esattezza del metodo, che passa da r = 2r + 1 a r = 2r − 1.
Sull’ambiente MATLAB® sono definite alcune funzioni utili al calcolo degli
integrali definiti:
I = integral(fun, a, b)
F = int(f,var)
78
incrementale:
f (x + h) − f (x)
lim ∀x ∈ [a, b]
h→0 h
Questo concetto è valido per tutti i metodi analitici per lo studio di funzione,
ma non è applicabile alle funzioni definite sui calcolatori, visto che essi sono affetti
da errori di risoluzione e men che meno hanno la capacità di astrarre il concetto
di limite. Per queste ragioni vanno trovati dei metodi di calcolo di tipo discreto,
in grado di approssimare le derivate.
Definizione 4.2.1 (Differenze finite in avanti). Usando una differenza finita po-
sitiva delle ascisse, si ottiene la differenza finita in avanti:
. f (x̄ + h) − f (x̄)
f ′ (x) ≃ δ+ f (x̄) = con h > 0
h
. f (x̄) − f (x̄ − h)
f ′ (x) ≃ δ− f (x̄) = con h > 0
h
Definizione 4.2.3 (Differenze finite centrate). Usando sia una differenza finita
positiva sia una differenza finita negativa delle ascisse, si ottiene la differenza finita
centrata:
. f (x̄ + h) − f (x̄ − h)
f ′ (x) ≃ δC f (x̄) = con h > 0
2h
Teorema 4.2.1 (Consistenza e convergenza metodi delle differenze finite). Tutti
e tre i metodi delle differenze finite sono consistenti e convergono alla derivata
prima:
lim δf (x̄) = f ′ (x̄)
h→0
79
Dimostrazione 4.2.1 (Consistenza e convergenza metodi delle differenze finite).
La dimostrazione del Teo. 4.2.1 è immediata: per verificare che il metodo sia
consistente (Def. 3 a pagina XII) si passa al limite delle differenze finite, che
coincide con la definizione di derivata prima: essendo il metodo consistente la
soluzione converge sempre.
Teorema 4.2.2 (Errori delle differenze finite in avanti e in indietro). Sia la fun-
zione f ∈ C 2 (Ix̄ ) gli errori delle differenze finite in avanti e all’indietro sono date
da:
h
E+ f (x̄) = f ′ (x̄) − δ+ f (x̄) = − f ′′ (ξ+ ) ξ+ ∈ [x̄, x̄ + h]
2
h
E− f (x̄) = f ′ (x̄) − δ− f (x̄) = − f ′′ (ξ− ) ξ− ∈ [x̄ − h, x̄]
2
Ciò implica che tutti polinomi di grado minore o uguale a uno vengono derivati
in modo esatto dalle differenze finite in avanti o all’indietro:
f (x̄ + h) − f (x̄)
E+ f (x̄) = f ′ (x̄) − δ+ f (x̄) = f ′ (x̄) −
h
f ′′ (ξ) 2
′ f (x̄) + f ′ (x̄)h + h − f (x̄) f ′′ (ξ) 2
E+ f (x̄) = f (x̄) − 2
=− h ξ ∈ [x̄, x̄ + h]
2 2
Teorema 4.2.3 (Errori delle differenze finite centrate). Sia la funzione f ∈ C 3 (Ix̄ )
gli errori delle differenze finite centrate sono dati da:
h [ ′′′ ]
EC f (x̄) = f ′ (x̄) − δC f (x̄) = − f (ξ+ ) + f ′′′ (ξ+ )
12
ξ+ ∈ [x̄, x̄ + h] ξ− ∈ [x̄ − h, x̄]
80
Dimostrazione 4.2.3 (Errori delle differenze finite centrate). L’errore di ap-
prossimazione è esprimibile come differenza tra la derivata nel punto e la sua
approssimazione:
EC f (x̄) = f ′ (x̄) − δC f (x̄)
1[ 1 1
EC f (x̄) = f ′ (x̄) − f (x̄) + f ′ (x̄)h + f ′′ (x̄)h2 + f ′′′ (ξ+ )h3 +
2h 2 6
1 h2 ] 1 [ ′′′ ]
− f (x̄) + f ′ (x̄)h − f ′′ (x̄)h2 + f ′′′ (ξ− )h3 = − f (ξ+ ) + f ′′′ (ξ− )
2 6 12
Per quanto riguarda gli ordini di convergenza di può facilmente affermare che
le differenze finite in avanti e all’indietro convergono linearmente, mentre quelle
centrate con ordine due.
h2 [ (4) ]
EC2 f (x̄) = f ′′ (x̄) − δC2 f (x̄) = −f (ξ+ ) + f (4) (ξ− )
24
ξ+ ∈ [x̄, x̄ + h] ξ− ∈ [x̄ − h, x̄]
81
Dimostrazione 4.2.4 (Differenze finite centrate per derivate seconde). Anco-
ra una volta l’errore è definito come differenza tra la derivata seconda e la sua
approssimazione numerica:
Questa volta si procedere con lo sviluppo in serie di Taylor per i termini f (x̄ + h)
e f (x̄ − h) fino al terzo ordine (è consentito farlo perché i termini di ordine due si
elidono a vicenda):
h2 [ (4) ]
EC2 f (x̄) = f ′′ (x̄) − δC2 f (x̄) = −f (ξ+ ) + f (4) (ξ− )
24
ξ+ ∈ [x̄, x̄ + h] ξ− ∈ [x̄ − h, x̄]
df = diff(f)
dove f è un vettore dove sono salvati tutte le valutazioni di una certa funzione,
mentre nel vettore df verranno calcolate tutte le differenze tra gli elementi suc-
cessivi del vettore di input. Questa funzione non calcola la differenza finita della
82
funzione: per ottenere l’algoritmo basta dividere ogni punto del vettore per la
distanza h che c’è fra ogni punto.
83
Capitolo 5
dy d2 y dn y
F (t, y, , 2 ,..., n ) = 0
dt dt dt
y (k) (t) t0
= y0,k ∀k = 0, 1, . . . , n − 1
84
funzione y(t) : I ⊂ R −→ R tale che:
y ′ (t) = f (t, y)
∀t ∈ I (5.1)
y(t ) = y
0 0
A questo punto si può integrare in modo definito usando come estremi di integra-
zione le condizioni iniziali del problema e le condizioni all’istante t:
∫ y(t) ∫ t
dy = g(τ )dτ ∀t ∈ I
y0 t0
85
5.1.1 Problema modello
Il problema modello è un particolare tipologia di ODE. di primo ordine, molto
facile da risolvere. Tale problema è molto spesso utilizzato per analizzare i metodi
numerici per la risoluzione di equazioni differenziali (da qui il nome). Esso è
definito come:
ẏ(t) = λy
∀t ∈ I ⊂ R λ ∈ R− (5.2)
y(t ) = y
0 0
y(t) = y0 eλ(t−t0 )
Dato che il numero reale λ è minore di zero, la soluzione dell’equazione y(t) tenderà
ad annullarsi per alti valori di t:
lim y(t) = 0
t→+∞
y(t) = y0 e[ℜ(λ)+iℑ(λ)](t−t0 )
[ ( ) ( )]
= y0 eℜ(λ)(t−t0 ) cos ℑ(λ)(t − t0 ) + i sin ℑ(λ)(t − t0 )
86
5.1.2 Soluzione all’equilibrio
Si consideri un problema di Cauchy del tipo:
f (t, y) = f (y) ∀t ∈ I
ȳ(t) f (ȳ(t)) = 0
ȳ1 = 0 ȳ2 = β
A questo punto sono distinguibili tre diversi intervalli: (−∞, ȳ1 ], (ȳ1 , ȳ2 ) e [ȳ2 , +∞).
87
divergere, poiché la derivata della funzione risulta monotona decrescente
ẏ = f (y) < 0.
• Allo stesso modo se il dato iniziale y0 = y(t0 ) > ȳ2 la funzione sarà an-
cora monotona decrescente, ma essendo la soluzione di equilibrio più pic-
cola della condizione iniziale, allora la soluzione y(t) tenderà a raggiungere
asintoticamente la soluzione ȳ2 .
Ossia se si ha un intervallo δ(ε) tale per cui ogni condizione iniziale y0 appartenente
all’intervallo non faccia si che la soluzione ad essa associata diverga. In formule,
questo concetto si esprime come:
lim y(t) = ȳ
t→+∞
88
5.1.3 Esistenza e unicità soluzione del problema di Cauchy
Il problema di Cauchy ammette l’esistenza di una sola soluzione solamente se
la funzione f (t, y) che descrive l’equazione differenziale ha certe proprietà. Si
consideri dunque il problema di Cauchy in forma normale (cioè quando la funzione
f esplicita la derivata di ordine massimo presente nell’equazione):
ẏ = f (t, y)
∀t ∈ I
t(t ) = y
0 0
∀y1 , y2 ∈ J
∃ L ∈ R+ |f (t, y1 ) − f (t, y2 )| < L|y1 − y2 |
∀t ∈ I
y1 →y2 y1 →y2
|y1 − y2 | −−−→ 0+ =⇒ |f (t, y1 ) − f (t, y2 )| −−−→ 0+ ∀t ∈ I
y1 →y2
f (t, y1 ) −−−→ f (t, y2 ) ∀t ∈ I
89
Teorema 5.1.4 (Esistenza e unicità locale per la soluzione del problema di Cau-
chy). Sia f una funzione f : I × J −→ R, se la funzione è continua su entrambi
gli argomenti e localmente Lipschitz nel secondo argomento allora la soluzione del
problema di Cauchy ad essa associata esiste ed è unica.
∀y1 , y2 ∈ R
∃ L ∈ R+ |f (t, y1 ) − f (t, y2 )| < L|y1 − y2 |
∀t ∈ I
∂f
|f (t, y1 ) − f (t, y2 )| = (t, ξ)(y1 − y2 ) ξ ∈ [y1 , y2 ]
∂y
∂f
|f (t, y1 ) − f (t, y2 )| ≤ max (t, y) |y1 − y2 | = L|y1 − y2 | ξ ∈ [y1 , y2 ]
t∈I ∂y
y∈R
90
è continua su tutto l’intervallo e limitata sull’intervallo chiuso e limitato [y1 , y2 ]
e quindi per il teorema di Weierstrass ammette massimo e minimo assoluti. Il
massimo è una costante di R e quindi si ritrova la definizione di funzione Lipschitz
(Def. 5.1.6).
in cui z(t) è la soluzione del problema perturbato, mentre δ0 e δ(t) sono le pertur-
bazioni sui dati.
91
allora: ∫t
φ(τ )dτ
w(t) ≤ g(t) e t0
Infine dato che la funzione C(t) è limitata sull’intervallo I = [t0 , tf ] si può definire
92
un massimo:
{ }
C = max C(t)
t∈I
F : I × Rm −→ Rm
f1 (t, y1 , y2 , . . . , ym )
f2 (t, y1 , y2 , . . . , ym )
F=
..
.
fm (t, y1 , y2 , . . . , ym )
93
A questo punto tutti i ragionamenti fatti per le equazioni differenziali di primo
ordine possono essere ripetuti tali e quali per le equazioni differenziali vettoriali.
Definizione 5.2.3 (Funzione vettoriale di Lipschitz). Si definisce funzione vetto-
riale localmente di Lipschitz una funzione tale che:
∑
m
y= Cj eλj (t−t0 ) vj
j=1
dove {λj }mj=1 e {vj }j=1 sono rispettivamente gli autovalori e gli autovettori asso-
m
ciati alla matrice A. Infine si ha che se tutti gli autovalori di A hanno parte reale
negativa, allora la soluzione è asintoticamente stabile e converge a zero:
{ }m
ℜ λj j=1 < 0 ⇐⇒ lim y = 0
t→+∞
94
sue derivate di primo e secondo ordine. In forma normale l’equazione si scrive:
f (t, y, ẏ) : I × R × R −→ R
v(t) = ẏ ∀t ∈ I
95
di ordine n. Si consideri ad esempio un problema di Cauchy del tipo:
y (n) = f (t, y, y ′ , . . . , y (n−1) )
y (t ) = yn−1
n−1 0
..
. ∀t ∈ I
y1 (t0 ) = y1
y (t ) = y
0 0 0
In questo caso è possibile riscrivere l’intero sistema tramite l’ausilio di unn vettore
di stato, del tipo:
( )T
z(t) = y (n−1) (t), y (n−2) (t), . . . , y ′ (t), y(t)
96
dove è una funzione incognita del tipo y(t) : I ⊆ Rm −→ Rm con m ≥ 1, le cui
componenti sono funzioni del tipo:
( )
y = x1 (t), x2 (t), . . . , xm (t)
97
equidistanti, dove h è il passo, mentre Nh il numero dei sottointervalli:
tf − t0
h=
Nh
tn = t0 + nh n = 0, 1, . . . , Nh
un − un−1
f (tn , un ) =
h
u0 = y0 (dato iniziale)
Per n = 0, 1, . . . , Nh − 1
un+1 = un + hf (tn , un )
98
Metodo di Eulero all’indietro
Questo metodo è del tutto analogo a quello predente, solo che utilizza le differenze
finite in indietro (Def. 4.2.2 a pagina 79) per approssimare la derivata prima.
( )
f tn , y(tn ) ≃ δ− f (tn )
u0 = y0 (dato iniziale)
Per n = 0, 1, . . . , Nh − 1
un+1 = un + hf (tn , un+1 )
Metodo di CrankNicolson
Questo metodo nasce come una media dei metodi di Eulero appena analizzati.
Questa volta la l’incremento delle funzione viene approssimato come:
∫ ∫
tn+1 un+1 ( )
∆y(tn+1 ) = ẏ(τ )dτ ≃ f τ, y(τ ) dτ = ∆un+1
tn un
h[ ]
un+1 − un = f (tn , un ) + f (tn+1 , un+1 )
2
99
Infine il metodo di Crank-Nicolson risulta:
u0 = y0 (dato iniziale)
Per n = 0, 1, . . . , Nh − 1
h[ ]
un+1 = un + f (tn , un ) + f (tn+1 , un+1 )
2
Anche in questo caso non si conosce il termine un+1 , per cui il metodo risulta
ancora una volta implicito.
Metodo di Heun
Finora si sono analizzati diversi metodi, alcuni più precisi (metodo Crank-Nicolson),
altri più veloci da risolvere (Eulero in avanti). Il metodo di Heun nasce con
l’obiettivo di rendere esplicito, e quindi più veloce, il metodo di Crank-Nicolson.
L’idea alla base è quella di stimare un+1 tramite il metodo di Eulero in avanti:
u∗n+1 = un + hf (tn , un )
h[ ]
un+1 − un = f (tn , un ) + f (tn+1 , u∗n+1 )
2
A questo punto l’espressione ottenuta non è più implicita, perché la funzione un+1
è stata approssimata a u∗n+1 allo step precedente.
L’algoritmo finale risulta dunque:
u0 = y0 (dato iniziale)
Per n = 0, 1, . . . , Nh − 1
u∗n+1 = un + hf (tn , un )
h[ ]
un+1 = un + f (tn , un ) + f (tn+1 , u∗n+1 )
2
100
5.3.2 Analisi degli errori e della stabilità
Come tutti i metodi numerici anche quelli costruiti per la risoluzione di ODE. sono
soggetti ad errori.
Si definisce errore alla discretizzazione di u la quantità:
ēn = y(tn ) − un
Metodo: un−1 −→ un
Calcolo errore: y(tn−1 ) −→ u∗n
In tal modo si riesce a distinguere l’errore dovuto alla consistenza del metodo,
creatosi al passo n, e l’errore imputabile alla stabilità del metodo.
y(tn ) − u∗n
τn (h) =
h
101
l’errore di troncamento globale si annulla al diminuire di h:
h→0
τ (h) −−→ 0+
τ (h) = O(hp )
h2
y(tn ) = y(tn−1 ) + hẏ(tn−1 ) + ÿ(ξn ) ξn ∈ [tn−1 , tn ]
2
h2 ( )
hτn (h) = y(tn−1 ) + hẏ(tn−1 ) + ÿ(ξn ) − y(tn−1 ) + hf (tn−1 , y(tn−1 ))
2
.
Tendendo presente che ẏ(tn−1 ) = f (tn−1 , y(tn−1 )) l’espressione si semplifica in:
h
τn (h) = ÿ(ξn )
2
h
τ (h) = D con: D = max |ÿ(ξn )|
2 n=1,...,Nh
102
u∗n − un :
Questo ragionamento può essere ripetuto n volte; sapendo inoltre che l’ultimo
termine della sommatoria è pari all’errore iniziale |ē0 | ed è nullo per le condizioni
iniziali del metodo si ottiene:
∑
n
|ēn | ≤ hτ (h) (1 + hL)k
k=0
∑
n
|ēn | ≤ hτ (h) (1 + hL)k ∼ hτ (h)(1 + hL)n ∼
k=0
enhL − 1 1 h
∼h τ (h) = (eL(tn −t0 ) − 1)D = C(D)h
hL L 2
103
Tabella 5.1: Nella tabella sono riportate le diverse tipologie, gli ordini di convergenza
e le ipotesi sulla funzione incognita per studiare gli errori di ogni metodo.
log(en1/en2 )
p=
log(hn1/hn2 )
Zero stabilità
z0 = y0 + δ0
Per n = 0, 1, . . . , Nh − 1
[ ]
zn+1 = zn + h f (tn , un ) + δn
104
Definizione 5.3.3 (Metodo Zero-stabile). Un metodo numerico per la risoluzione
di equazioni differenziali ordinarie si dice zero-stabile nell’intervallo I = [t0 , tf ]
limitato se esistono h0 , C, ε0 ∈ R+ tali che:
C = eL(tf −t0 ) − 1
tf →+∞
C −−−−→ +∞
Assoluta stabilità
In altre parole si può dire che l’assoluta stabilità sia una generalizzazione della
zero stabilità a domini molto più ampi.
105
che:
lim un = 0
n→+∞
h ∈ (0, hmax ]
Nei prossimi paragrafi viene analizzata l’assoluta stabilità dei quattro metodi
presentati finora.
=⇒ |1 + hλ| < 1
−1 < 1 + hλ < 1
−2 < hλ < 0
106
Quindi affinché il metodo du Eulero avanti sia assolutamente stabile deve essere
rispettata la seguente condizione:
( )
2
h∈ 0, hmax =
|λ|
Passando al limite:
y0
lim un = lim =0
n→+∞ n→+∞ (1 + hλ)n
=⇒ |1 − hλ| > 1
h[ ]
un+1 = un + f (tn , un ) + f (tn+1 , un+1 )
2
h
un+1 = un + λ[un + un+1 ]
( ) 2 ( )
h h
1 − λ un+1 = 1 + λ un
2 2
( )n+1
1 + h2 λ 1 + h2 λ
un+1 = un = y0
1 − h2 λ 1 − h2 λ
107
Passando al limite:
( )n
1 + h2 λ
lim un = lim y0 = 0
n→+∞ n→+∞ 1 − h2 λ
1 + h2 λ
=⇒ <1
1 − h2 λ
1 + h2 λ 2 + hλ 2 − h|λ|
= = <1 ∀h ∈ R+
1 − h2 λ 2 − hλ 2 + h|λ|
Passando al limite:
1
lim un = lim (2 + 2hλ + (hλ)2 )n+1 y0 = 0
n→+∞ n→+∞ 2n+1
2 + 2hλ + (hλ)2
=⇒ <1
2
−2 < 2 + 2hλ + (hλ)2 < 2
−4 < 2hλ + (hλ)2 < 0
108
Tabella 5.2: Funzioni di stabilità (a variabile complessa z = hλ) dei quattro metodi
analizzati.
Heun RH (z) = 1 + z + 12 z 2
R : C −→ C
[ ]n
un = R(hλ y0
|R(hλ)| < 1
A titolo di esempio, in Fig. 5.1 sono riportate le aree di assoluta stabilità dei
109
metodi analizzati fino ad adesso.
∂f
(t, y) < 0 ∀t ∈ [t0 , +∞)
∂y
2 ∂f
0 < h < hmax = λmax = max (t, y)
λmax t∈[t0 , +∞) ∂y
u0 = y0 (Dato iniziale)
Per n = 0, 1, 2, . . . , Nh − 1
∑s
un+1 = un + h bi ki
i=1
110
2 2
ℑ(z) ℑ(z)
1 1
ℜ(z) ℜ(z)
−3 −2 −1 1 −1 1 2 3
−1 −1
−2 −2
(a) Il metodo di Eulero avanti risulta as- (b) Il metodo di Eulero indietro risulta
solutamente stabile per tutti i valori di assolutamente stabile per tutti i valori di
z = hλ che appartengono all’area di un z = hλ che non appartengono all’area di
cerchio di raggio uno centrato in z0 = −1. un cerchio di raggio uno centrato in z0 = 1.
2 2
ℑ(z) ℑ(z)
1 1
ℜ(z) ℜ(z)
−3 −2 −1 1 −3 −2 −1 1
−1 −1
−2 −2
(c) Il metodo di Crank-Nicolson risulta (d) Il metodo di Heun risulta assolutamen-
assolutamente stabile per tutti i valori di te stabile per tutti i valori di z = hλ che
z = hλ che appartengono al semipiano appartengono all’area un’ellisse centrata in
complesso avente parte reale negativa. z0 = −1.
Figura 5.1: In figura sono rappresentate le quattro regioni di assoluta stabilità per i
metodi analizzati finora.
111
Si può notare facilmente che la funzione f non compare nell’algoritmo, almeno in
maniera esplicita, proprio perché le valutazioni sono fatte su s stadi che non hanno
uguale misura tra loro. Analizzando più a fondo l’espressione si può notare che la
sommatoria presente ad ogni passaggio può essere vista come il prodotto scalare di
due vettori k, b ∈ Rs , in cui il primo contiene le valutazioni della funzione f (t, y)
in opportuni punti, mentre il secondo i pesi di ognuna di esse. Nello specifico il
vettore k è calcolato come:
( ∑
s )
ki = f tn + ci h, un + aij kj ∀i = 1, 2, . . . , s
j=1
c A
bT
0 0 0
1 1 0
1/2 1/2
112
Calcolo dei coefficienti ki :
h
un+1 = un + [fn + f (tn , un + fn )]
2
s2 < s 1 ⇐⇒ usn+1
1
− usn+1
2
> ε = toll
In caso contrario si può optare per hn+1 > hn , per risparmiare tempo di calcolo.
In MATLAB® sono implementate diverse funzioni atte a risolvere le equazioni
differenziali ordinarie; in particolare si segnalano:
[t, y] = ode23(fun, T, y0)
[t, y] = ode45(fun, T, y0)
Entrambe le funzioni restituiscono i vettori t e y che contengono i tempi e le
valutazioni della risoluzione approssimata e richiedono in ingresso la funzione fun,
l’intervallo di integrazione T e la condizione iniziale y0. Ovviamente le due funzioni
non sono equivalenti, ma utilizzano due diversi ordini di metodo Runge-Kutta
espliciti: la prima ha un tempo di esecuzione molto minore, ma meno precisa dal
punto di vista degli errori.
113
5.3.4 Theta-Metodi
Anche in questo caso il punto di partenza è la discretizzazione del dominio tem-
porale I in un numero Nh di punti equidistanti.
Questa volta si introduce un parametro ϑ ∈ [0, 1] la cui funzione è quella di
pesare le valutazioni della funzione f tra i due passi consecutivi tn e tn+1 . I metodi
risultano quindi:
u0 = y0 (dato iniziale)
Per n = 0, 1, . . . , Nh − 1
[ ]
un+1 = un + h (1 − ϑ)f (tn , un ) + ϑf (tn+1 , un+1 )
1. (ϑ = 0) Eulero avanti:
un+1 = un + hf (tn , un )
2. (ϑ = 1) Eulero indietro:
h[ ]
un+1 = un + f (tn , un ) + f (tn+1 , un+1 )
2
114
5.3.5 Risoluzione di sitemi differenziali ordinari
Tutti i metodi analizzati sono applicabili alla risoluzione di sistemi di funzioni
differenziali ordinarie.
Per prima cosa si riduce il problema a un problema di Cauchy vettoriale del
primo ordine, come visto in precedenza; poi si applica uno dei metodi visti prima.
Ad esempio i ϑ-metodi diventano:
u0 = y0 (dato iniziale)
Per n = 0, 1, . . . , Nh − 1
[ ]
un+1 = un + h (1 − ϑ)F(tn , un ) + ϑF(tn+1 , un+1 )
Nel particolare caso dei ϑ-metodi è possibile fare una riflessione in merito alla
funzione F(t, y); supponiamo che per ipotesi sia possibile scrivere la funzione come:
F(tn , y) = Ay + g(t)
Questa espressione è di fatto un sistema lineare: tutto quello che c’è a destra
dell’uguale è costante e costituisce il termine noto, mentre (I − hϑA) è la matrice
del sistema e un+1 rappresenta l’unica incognita dell’espressione.
115
Teorema 5.3.2 (Effetto autovalori matrice di stato). Se gli autovalori della ma-
trice di stato A hanno tutti parte reale negativa allora:
{ }m
ℜ λi (A) <0 =⇒ lim y(t) = 0
t→+∞
i=1
Come fatto precedentemente questo criterio può essere esteso a tutte le funzioni
F(t, y) generiche, andando a considerare gli autovalori della matrice Jacobiana. In
particolare, il metodo risulta stabile se:
[ ] ∂fi
JF (t) ij = (t, y)
∂yj
( )
λi JF (t) ∈ A ∀i = 1, 2, . . . , m ∀t ∈ I
Metodo Leap-Frog
116
Sia dunque un l’approssimazione della funzione al passo n e vn l’approssimazione
della sua derivata prima:
h2
un+1 = un + hvn + f (tn , un , vn )
2
1
vn = (un+1 − un−1 )
2h
1
vn+1 = (un+2 − un )
2h
A questo punto si sostituisce un+2 , due volte, con il valore appena calcolato:
( )
1 h2
vn+1 = un+1 + hvn+1 + f (tn+1 , un+1 , vn+1 ) − un
2h 2
( 2
)
1 h h2
= un + hvn + f (tn , un , vn ) + hvn+1 + f (tn+1 , un+1 , vn+1 ) − un
2h 2 2
117
Infine si ricava vn+1 :
h( )
vn+1 = vn + f (tn , un , vn ) + f (tn+1 , un+1 , vn+1 )
2
u0 = y0 (Dato iniziale)
v0 = y1 (Dato iniziale)
Per n = 0, 1, . . . , Nh − 1
h2
un+1 = un + hvn + f (tn , un , vn )
2
h( )
vn+1 = vn + f (tn , un , vn ) + f (tn+1 , un+1 , vn+1 )
2
Dall’algoritmo si capisce subito il perché del nome del metodo (letteralmente signi-
fica balzo della rana): ad ogni passaggio calcola la funzione nel punto e la derivata
nel nodo immediatamente successivo; successivamente torna indietro a calcolare la
funzione. . .
Si nota subito che la prima equazione è lineare in un , quindi non ci sono partico-
lari problemi; purtroppo non è così per la seconda: il termine vn+1 non è necessaria-
mente lineare (dipende dalla funzione f ) quindi potrebbe richiedere l’applicazione
di uno dei metodi di risoluzione di equazioni non lineari.
In ogni caso, il metodo Leap-Frog ha ordine di convergenza p = 2 ed è sempre
assolutamente e incondizionatamente stabile.
118
Metodo di Newmark
Il metodo di Newmark nasce come estensione dei ϑ-metodi ai problemi del secondo
ordine. Assegnati ξ, ϑ ≥ 0 si definisce l’algoritmo come:
u0 = y0 (Dato iniziale)
v0 = y1 (Dato iniziale)
Per n = 0, 1, . . . , Nh − 1
h2 [ ]
un+1 = un + hvn + 2ξfn+1 + (1 − ξ)fn
[ 2 ]
vn+1 = vn + h ϑfn+1 + (1 − ϑ)fn
• il metodo è esplicito se ξ = ϑ = 0;
1
• il metodo coincide con Leap-Frog se ϑ = 2
e ξ = 0.
1
2ξ ≥ ϑ ≥
2
119
Parte II
120
Capitolo 6
Introduzione
Una delle distinzioni più importanti che si possono fare per analizzare le equa-
zioni differenziali è quella tra equazioni lineari, in cui la funzione F è lineare
121
rispetto alla funzione u e a tutte le sue derivate, e equazioni non lineari. Per
motivi di semplicità verranno proposte solamente problemi di tipo lineare.
• l’equazione differenziale;
• una o più funzioni del tempo che rappresentino le condizioni al contorno per
la variabile spaziale.
122
Se invece la funzione incognita è definita su una variabile spaziale di tipo vettoriale
allora la condizione diventa:
Nel caso in cui la funzione g(t) (o più di una nel caso a più dimensioni)
sia costantemente nulla, si dice che le condizioni al contorno di Dirichlet sono
omogenee.
Allo stesso modo se la funzione vettoriale dipende da una variabile spaziale di tipo
vettoriale allora la condizione di Neumann diventa:
123
infinito) e la funzione u(x, t); nel caso a una sola dimensione questo diventa:
∂u (0, t) = γ (U − u(0, t))
∂x
∀t ∈ (t0 , tf ) γ ∈ R
(L, t) = γ (U − u(0, t))
∂u
∂x
Nel caso il problema sia definito per funzioni spaziali a dominio vettoriale la
condizione diventa:
( )
∇u(x, t) · n̂ = γ U − u(x, t) ∀t ∈ (t0 , tf ) γ ∈ R
∂Ω
∃ T ∈ R+ f (x) = f (x + T ) ∀x ∈ R
124
Affinché si possa definire lo sviluppo in serie di Fourier di una certa funzione
f è necessario che questa sia dotata di una minima regolarità. La condizione
sufficiente a soddisfare l’ipotesi e che la funzione appartenga allo spazio L2 .
a0 ∑ ( )
+∞
f (x) = + an cos(nx) + bn sin(nx)
2 n=1
∑
+∞ ∫ T
1
f (x) = ck e ink
cn = f (x)e−inx dx
n=−∞
2T −T
1 1 1
c 0 = a0 cn = (an − ibn ) c−n = c̄n = (an + ibn )
2 2 2
a0 = 2c0 an = cn + c−n bn = i(cn − c−n )
125
Teorema 6.3.2 (Proprietà funzioni in L2 ). Una funzione che appartiene allo
spazio L2 su un intervallo (0, L) dispone delle seguenti proprietà:
2. il limite:
∫ L [ ]2
lim ∥fN (x) − f (x)∥ = lim fN (x) − f (x) dx = 0
N →+∞ N →+∞ 0
4. I limiti:
lim an = lim bn = lim |cn | = 0
n→+∞ n→+∞ n→+∞
126
Capitolo 7
Problemi di diffusione
∂u
(x, t) − D∇2 u(x, t) = f (x, t)
∂t
k=1
∂k ∂k
127
7.2 L’equazione del calore
Si consideri ad esempio di voler studiare il comportamento termico di un filo me-
tallico, avente una dimensione pari a L molto maggiore delle altre due (problema
monodimensionale). Si immerga in contatto tale filo in un ambiente a temperatura
nota e costante, pari a T0 , per un tempo sufficientemente lungo affinché si instauri
l’equilibrio termico:
u(x) = T0 ∀x ∈ [0, L]
Si metta a contatto il filo con due corpi aventi temperatura pari a T1 e T2 : dopo
un po’ di tempo si osserva che il filo tende ad assumere una distribuzione di tipo
lineare, i cui massimi e minimi sono posti in corrispondenza degli estremi. Il
caso appena descritto è un tipico comportamento che l’equazione di diffusione può
analizzare e prevedere senza particolari difficoltà.
Teorema 7.2.1 (della divergenza). Si consideri un volume V ⊆ Rn compatto e
delimitato da una superficie di frontiera liscia ∂S. Se la funzione F è una funzione
differenziabile con continuità su tutto il volume, si ha che:
∫ ∫
∇ · FdΩ = F · n̂dσ
V ∂V
dove EV è l’energia totale del corpo, mentre ρ la sua densità. In tal modo il tasso
di variazione è pari a:
∫ ∫
d
ĖV = ρe(x, t)dΩ = ρet (x, t)dΩ
dt V V
128
volume ∂V . Sotto tali ipotesi l’energia entrante (o uscente se di segno negativo) è
pari a: ∫ ∫
Ein = − q · n̂dσ = − ∇ · qdΩ
∂V V
Dato che non si è fatta alcun genere di ipotesi sul dominio V l’uguaglianza risulta
verificata anche nel caso in cui si rimuova l’integrazione:
129
Problema di Cauchy-Dirichlet Si definisce il problema di Cauchy-Dirichlet il
problema differenziale alle derivate parziali del tipo:
(x, t) ∈ (0, L) × (0, tf )
u = Duxx
t
u(0, t) = g1 (t) u(L, t) = g2 (t) t ∈ (0, tf )
u(x, 0) = ϕ(x) x ∈ (0, L)
. { }
QT = (x, t) ∈ (0, L) × (0, tf )
anche se la funzione che descrive la condizione iniziale di Cauchy ϕ(x) ̸∈ C 0 ((0, L)).
Tra i quattro problemi ce n’è uno che ha una proprietà in più degli altri: quello
di Neumann. Esso infatti, quando ha le condizioni al contorno omogenee, oltre ad
avere una sola soluzione, ha valor medio integrale sullo spazio costante nel tempo.
130
Dimostrazione 7.3.1 (Valor medio problema Cauchy-Neumann omogeneo). An-
dando a moltiplicare l’equazione differenziale omogenea per una funzione generica
v(x) ̸= 0:
∂u ∂ 2u
v = Dv 2
∂t ∂x
Integrando su tutto il dominio spaziale si ha che:
∫ L ∫ L
∂u ∂ 2u
v dx = D v dx
0 ∂t 0 ∂x2
Dato che l’espressione è valida per tutte le funzioni v(x) ̸= 0 si può scegliere una
funzione costante, ad esempio v(x) ≡ 1, ∀x ∈ (0, L); in tal modo si annulla il
secondo termine dell’equazione:
∫ L ∫ L
∂u d
dx = u(x, t)dx = 0
0 ∂t dx 0
131
Sostituendo la definizione nel passaggio precedente si ottiene:
d
L um = 0
dx
Da cui segue che il valore medio sullo spazio della soluzione u(x, t) è una funzione
tempo invariante.
Si supponga che per assurdo esistano due funzioni diverse u e v che siano
soluzioni del problema differenziale. Essendo l’equazione di diffusione lineare, una
qualsiasi combinazione lineare delle soluzioni risulta essere ancora soluzione del
problema.
∃ w(x, t) = u(x, t) − v(x, t) ̸= 0 ∀(x, t) ∈ QT
132
Sostituendo l’espressione nel problema, si ottiene il seguente sistema:
(x, t) ∈ QT
w (x, t) = Dwxx (x, t)
t
w(x, 0) = ϕ(x) − ϕ(x) ≡ 0 x ∈ (0, L)
Condizione omogenea di D, N, R, M.
∂w ∂ 2w
(x, t) = D 2 (x, t)
∂t ∂x
Cioè è pari alla derivata dell’energia della funzione w(x, t) cioè alla derivata dell’in-
tegrale di una funzione al quadrato, cioè a quantità sempre positiva o al massimo
nulla. Al contrario il secondo integrale può essere riarrangiato sfruttando la regola
di differenziazione del prodotto di funzioni (integrazione per parti):
133
∫ L
∂ 2w
II = D w
dx
∂x2
0
∫ L ( ) ∫ L
PP ∂ ∂w ∂w ∂w
=D w dx − D dx
0 ∂x ∂x 0 ∂x ∂x
[ ]L ∫ L [ ]2
∂w ∂w
=D w −D dx
∂x 0 0 ∂x
Per il teorema del confronto semplice la derivata dell’energia è sempre nulla, quindi
la funzione w(x, t), sarà una costante nel tempo.
Si consideri ora, l’energia che si ha all’istante t = 0, cioè quando la funzione
deve essere pari alla condizione iniziale del problema:
∫
1 L ( )2
E(0) = w(x, 0) dx = 0
2 0
Da questo segue che, fra le infinite costanti possibili l’unica accettabile è quella
nulla:
E(t) ≡ 0 ∀t ∈ (0, tf )
Il valor medio dell’energia della funzione w(x, t) sullo spazio è a sua volta nullo:
∫ L ( )2
w(x, t) dx ≡ 0 ⇐⇒ w(x, t) ≡ 0 ∀(x, t) ∈ QT
0
134
E questo a sua volta è verificato solo se:
Cioè le due funzioni, diverse per ipotesi, devono essere uguali puntualmente, che
è un assurdo. La soluzione a tutti i problemi diffusione con condizioni di Cauchy
e bordo, ammesso che esista, risulta unica.
QT = Ω × (t0 , tf )
135
Dimostrazione 7.4.2 (Unicità soluzione multidimensionale). Sia dato un pro-
blema di diffusione omogeneo su uno spazio multidimensionale, avente una fra le
quattro tipologie di condizione al contorno analizzate.
u (x, t) = D∇2 (x, t) (x, t) ∈ QT
t
u(x, 0) = ϕ(x) x∈Ω
Condizione di D, N, R, M.
Come prima si suppone che per assurdo esistono due soluzioni distinte u(x, t) e
v(x, t) dello stello problema di diffusione. Essendo l’equazione di diffusione lineare
ogni combinazione lineare delle due soluzioni risulta deve per forza essere anch’essa
soluzione all’equazione. Si consideri quindi la soluzione w(x, t) definita come:
Con tale funzione le condizioni al bordo risultano omogenee, qualsiasi sia la loro
tipologia.
Con passaggi analoghi al caso monodimensionale (moltiplicazione per la fun-
zione w e integrazione su tutto il dominio spaziale) si arriva a:
∫ ∫
wwt dω = D w∇2 wdω
Ω Ω
Per via delle condizioni al contorno omogenee il termine di flusso è sempre nullo:
nello specifico Dirichlet annulla la funzione, Neumann e Robin la sua derivata
136
(rappresentata dal gradiente). Da questo deriva che:
∫
II = − ∥∇w∥2 dω ≤ 0 ∀t ≥ 0
Ω
e questo rappresenta un assurdo. Si conclude quindi che non esistono due soluzioni
distinte all’equazione di diffusione.
137
Teorema 7.5.1 (Proprietà serie di Weierstrass). Data una serie di funzioni di
∑
Weierstrass n fn (x) su un intervallo I, in cui tutte le funzioni fn (x) sono continue
sull’intervallo allora:
∑
+∞
f (x) = fn (x) ∈ C 0 (I)
n=1
• il limite delle serie esiste finito ed è uguale alla serie dei limiti:
∃ lim fn (x) = ln ∈ R
x→x0
∑
+∞
∃ lim f (x) = ln = l ∈ R
x→x0
n=1
∑
+∞
′
f (x) = fn′ (x)
n=1
Ora si dispone di tutti gli strumenti necessari per la risoluzione del problema
di diffusione. Per farlo si adotta la tecnica della separazione delle variabili, che
prevede che la soluzione al problema sia pari al prodotto di due funzioni distinte,
ognuna delle quali è una funzione di una sola variabile:
138
• sono lineari;
• sono omogenee;
ut − Duxx = 0
∂U ∂ 2U
=D 2
∂t ∂x
∂(XT ) ∂ 2 (XT )
=D
∂t ∂x2
XT = DT X ′′
′
1 T′ X ′′
=
DT X
139
Naturalmente questa espressione deve valere in ogni punto del dominio QT , poiché
la stessa equazione è definita su esso. Affinché tale relazione sia soddisfatta è
necessario che entrambi i rapporti siano pari ad una costante, per via del fatto che
le funzioni T e X non hanno nessuna variabile in comune.
1 T′ X ′′
= =k∈R
DT X
T (t) = eDkt
Il problema è ora capire che valore ha la costante k, per farlo si utilizza il secondo
problema differenziale, cioè quello relativo alla funzione X(x):
X ′′ (x) − kX(x) = 0
in funzione del valore che assume la costante si hanno dei comportamenti diversi
della soluzione; nello specifico si distinguono tre casi:
X ′′ (x) − λ2 X(x) = 0
µ2 − λ2 = 0 =⇒ µ = ±λ
140
In tal caso la soluzione è pari a:
1 1
−λL +λL
= e+λL − e−λL ̸= 0 ∀λ > 0
e e
Quindi per il teorema di Cramer esiste una sola soluzione; essendo inoltre il
sistema omogeneo si ha che la soluzione è data da:
c1 = c2 = 0
X(x) ≡ 0 x ∈ (0, L)
2. La costante è nulla: k = 0
Questo semplifica notevolmente il problema differenziale, che si riduce a:
X ′′ (x) = 0
X(x) = c1 x + c2
141
Come prima si impongono le condizioni al bordo, ottenendo il sistema lineare:
[ ]{ }
X(0) = c = 0 0 1 c1
2
=⇒ =0
X(L) = c L + c = 0 L 1 c2
1 2
Ancora una volta il determinante della matrice non può mai essere nullo:
0 1
= −L ̸= 0
L 1
c1 = c2 = 0
X ′′ (x) + λ2 X(x) = 0
µ2 + λ2 = 0 =⇒ µ = ±iλ
In tal caso la soluzione è pari alla somma di due esponenziali complessi, che
per l’identità di Eulero si possono riscrivere come:
142
lineare:
[ ]{ }
X(0) = c = 0 1 0 c1
1
=⇒ =0
X(L) = c cos(λL) + c sin(λL) = 0 cos(λL) sin(λL) c2
1 2
In questo caso il determinante della matrice può essere nullo per alcuni
particolari valori di λ:
1 0 nπ
= sin(λL) = 0 ⇐⇒ λ= n∈N
cos(λL) sin(λL) L
Un (x, 0) = ϕ(x)
∑
+∞ ( )
nπ
ϕ(x) = bn sin x
n=1
L
143
Da cui segue, per il principio di identità, che tutti i coefficienti cn , ancora incogniti,
sono pari ai coefficienti di Fourier bn e quindi la funzione U (x, t) è pari a:
∑
+∞ ( )
−D( nπ )2 t nπ
U (x, t) = bn e L sin x
n=1
L
Adesso però è necessario verificare che la soluzione U (x, t) appena costruita con-
verga alla soluzione del problema u(x, t) e che soddisfi l’equazione differenziale.
Per risolvere il primo problema, si deve maggiorare la serie, in modo da lavorare
con termini più semplici da trattare:
( )
−D( nπ nπ
≤ |bn |e−D( L )
nπ 2
)2 t t
bn e L sin x
L
q<1
Tale ipotesi è senz’altro verificata, per cui la serie U (x, t) converge alla funzione
u(x, t) per tutti i punti appartenenti al dominio QT .
L’ultimo problema da affrontare consiste nel verificare la soluzione u(x, t) soddi-
sfa l’equazione differenziale. Sappiamo che sicuramente la funzione U (x, t) verifica
l’equazione, perché è stata costruita a partire da essa, e tende alla funzione u(x, t).
Rimane solo da vedere se anche la funzione u(x, t) soddisfa l’equazione di diffu-
sione. In altre parole bisogna verificare che la seguente espressione risulti essere
un’identità:
( )∑
+∞
∗ ?
∑
+∞
( ) ∗
∂t − D∂xx
2
bn Un = bn ∂t − D∂xx 2
Un
n=1 n=1
in cui il termine Un∗ è pari alla funzione Un divisa per il suo coefficiente bn .
Poco fa si è dimostrato che tutti i termini Un sono più piccoli di una specifica
144
costante qn , per cui collezionando queste costanti, si può determinare una succes-
∑
sione geometrica convergente: per questo motivo la serie n bn Un∗ è una serie di
Weierstrass (Def. 7.5.1) e l’uguaglianza vista prima risulta sempre verificata.
La soluzione al problema di diffusione omogeneo, con condizioni al contorno di
Dirichlet omogenee, definito dall’Eq. (7.2), è dunque:
∑
+∞ ( )
−D( nπ )2 t nπ
u(x, t) = bn e L sin x
n=1
L
dove i coefficienti bn sono i coefficienti della serie di Fourier della funzione ϕ(x):
∫ L ( )
2 nπ
bn = ϕ(x) sin x dx
L 0 L
lim u(x, t) = 0
t→+∞
145
sicuramente derivata temporale nulla, poiché non dipende dal tempo:
∂us
=0 ∀x ∈ (0, L) (7.4)
∂t
∂ 2 us
=0 (7.5)
∂x2
∂ ∂2
(us (x) + ū(x, t)) = D 2 (us (x) + ū(x, t))
∂t ∂x
∂ ∂2
ū(x, t) = D 2 ū(x, t)
∂t ∂x
= a − ū(0, t)
146
Quindi dato che il problema impone che u(0, t) = a, si ha che:
ū(0, t) = a − a = 0
∑
+∞ ( )
−D( nπ )2 t nπ
u(x, t) = b̄n e L sin x + us (x)
n=1
L
dove i coefficienti b̄n sono i coefficienti della serie di Fourier della funzione ϕ̄(x):
∫ L ( )
2 nπ
b̄n = ϕ̄(x) sin x dx
L 0 L
147
quindi:
u − Duxx = f (x) (x, t) ∈ (0, L) × (0, tf )
t
u(0, t) = u(L, t) = 0 t ∈ (0, tf )
u(x, 0) = ϕ(x) x ∈ (0, L)
Si procede ancora con una soluzione di tentativo, del tipo:
∑
+∞ ( )
nπ
u(x, t) = cn (t) sin x
n=1
L
∑
+∞ ) (
nπ
ut (x, t) = c′n (t) sin
x
n=1
L
( )
n2 π 2 ∑
+∞
nπ
uxx (x, t) = − 2 cn (t) sin x
L n=1 L
+∞ [
∑ ] ( )
′ n2 π 2 nπ
cn (t) − 2 cn (t) sin x = f (x)
n=1
L L
∑
+∞ ( ) ∫ L ( )
nπ 2 nπ
f (x) = kn sin x kn = f (x) sin x dx
n=1
L L 0 L
n2 π 2
c′n (t) − cn (t) = kn n ∈ N∗
L2
148
Il lato positivo di tutte queste equazioni differenziali è che sono tutte funzioni della
sola variabile t. La generica soluzione è data da:
n2 π 2 2
t L2
cn (t) = An e L2 + kn (7.6)
n2 π 2
Infine, come già fatto per i precedenti problemi, si sviluppa in serie di Fourier la
funzione ϕ(x) e la si pone uguale alla soluzione trovata, in modo da calcolare tutti
i coefficienti An ancora incogniti:
∑
+∞ ( ) ∑
+∞ ( )
nπ nπ
u(x, 0) = cn (0) sin x = bn sin x
n=1
L n=1
L
L2
An = bn − kn ∀n ∈ N∗
n2 π 2
+∞ [
∑ ( )] ( )
n2 π 2 2 L2 n2 π 2 2 nπ
u(x, t) = bn e L2
t
+ 2 2 kn 1 − e L 2 t
sin x
n=1
nπ L
149
Condizioni al contorno omogenee
Come sempre si parte dal facile: funzione omogenea con condizioni al contorno
omogenee. Il problema differenziale risulta quindi essere:
u − Duxx = 0 (x, t) ∈ (0, L) × (0, tf )
t
ux (0, t) = ux (L, t) = 0 t ∈ (0, tf )
u(x, 0) = ϕ(x) x ∈ (0, L)
1 T′ X ′′
= =k∈R
DT X
X ′′ (x) = λ2 X(x)
150
per cui esiste un’unica soluzione:
c1 = c2 = 0
X ′′ (x) = 0
X(x) = c1 x + c2
151
Sostituendo le condizioni al contorno di Neumann si trova:
X ′ (0) = c λ = 0
1
X ′ (L) = c λ cos(λx) − c λ sin(λL) = 0
1 2
nπ
sin(λx) = 0 ⇐⇒ Lλ = nπ ⇐⇒ λ= n ∈ N∗
L
∑
+∞ ( )
nπ
X(x) = cn cos x
n=1
L
nπ
T ′ (t) = Dλ2 T (t) = D T (t)
L
nπ
T (t) = eD L t
152
rimane solo da applicare la condizione iniziale di Cauchy, data dalla funzione ϕ(x)
e da verificare che la soluzione totale U (x, t) trovata soddisfi ancora l’equazione di
diffusione. Espandendo la funzione ϕ(x) tramite serie di Fourier di soli coseni, si
ottiene:
∑ ( ) ( )
a0 ∑
+∞ +∞
nπ nπ
U (x, 0) = co + cn cos x = + an cos x = ϕ(x)
n=1
L 2 n=1
L
a0
lim u(x, t) =
t→+∞ 2
Tale concetto è in pieno accordo con tutto quello precedentemente dimostrato nella
Sez. 7.3 a pagina 129.
153
Condizioni al contorno costanti
Anche in questo caso non si può risolvere direttamente il problema con il meto-
do della separazione delle variabili, ma bisogna per prima cosa rendere omogenee le
condizioni al bordo. Si consideri la soluzione al problema come una somma di due
funzioni: una funzione che descrive il transitorio ū(x, t), che quindi tende ad annul-
larsi nel tempo, e una soluzione stazionaria us (x) che descrive il comportamento
del sistema a tempo infinito.
∂us ∂ 3 us
−D 2 =0
∂t ∂x
La soluzione us è una soluzione che per ipotesi non dipende dal tempo, per cui
l’equazione di diffusione si riduce a:
∂ 3 us
=0
∂x2
154
in (0, L), ma una parabola del tipo:
us (x) = c1 x2 + c2 x + c3
b−a 2
us (x) = x + ax
2L
∂ 3 us b−a
2
=
∂x L
b−a 2
us (x, t) = x + ax + ct t
2L
155
rimane dalla derivata seconda spaziale sia annullata dalla costante ct :
∂us ∂ 2 us
−D 2 =0
∂t ∂x
b−a
ct − D =0
L
b−a
=⇒ ct = D
L
b−a 2 b−a
us (x, t) = x + ax + D t
2L L
Trovata la prima soluzione si passa allo studio del vero e proprio transito-
rio. Dalla relazione data dall’Eq. (7.8) e sostituendola nel problema descritto
dall’Eq. (7.7) si ottiene:
∂( ) ∂2 ( )
us (x, t) + ū(x, t) − D 2 us (x, t) + ū(x, t) = 0
∂t ∂x
∂ ū ∂ 2 ū
−D 2 =0 ∀(x, t) ∈ QT
∂t ∂x
156
con condizioni di Neumann omogenee:
∂ ū 2
− D ∂∂xū2 = 0 ∀(x, t) ∈ QT
∂t
Teorema 7.5.2 (Principio del massimo e del minimo). Sia u(x, t) una funzio-
ne continua su tutto il dominio e differenziabile due volte rispetto alla variabile
spaziale x e una su quella temporale t, in tutti i punti interni.
157
• ut − ∇2 u = q ≤ 0 in QT allora il massimo della soluzione si trova sulla
frontiera parabolica:
{ } { }
max u(x, t) = max u(x, t)
(x,t)∈Q̄T (x,t)∈∂P QT
wt − D∇2 w = q − ε < 0
∇2 w(x0 , t∗ ) ≤ 0
158
Inoltre la derivata temporale è:
w (x , t∗ ) = 0 se: t∗ < T − ε
t 0
w (x , t∗ ) ≥ 0 se: t∗ = T − ε
t 0
Per passare dal massimo della funzione sulla sotto regione Q̄T −ε al massimo della
funzione sull’intera regione Q̄T è sufficiente far tendere ε → 0:
max u −→ max u
Q̄T −ε Q̄T
ε→0
max u + T ε −−→ max u
∂P Q T ∂P Q T
Partendo dal principio del massimo si possono fare alcune importanti conside-
razioni sull’equazione del calore. Per prima cosa si considerino due funzioni diverse
u, w ∈ C 2,1 (QT ) ∩ C 0 (Q̄T ), che soddisfano i problemi:
ut − D∇2 u = f1 wt − D∇2 w = f2
159
Inoltre vale:
Questo permette di fare dei confronti tra le diverse funzioni, anche senza sapere le
loro espressioni analitiche. Inoltre è anche possibile affermare che il problema di
diffusione di Cauchy-Dirichlet dipende in modo continuativo dai dati:
u(x, t) ∂Ω
= g(t) w(x, t) ∂Ω
= g(t) + ε f2 = f1 + ε
Di fatto si sta dicendo che w è uno stato perturbato di u rispetto ai dati al bordo.
Applicando le Eq. (7.10) e Eq. (7.11):
max |u − w| ≤ ε + T ε = (1 + T )ε
(x,t)∈Q̄T
Questo indica che le perturbazioni dei dati in un generico istante sono limitate da
quelle che si hanno all’istante iniziale e dal tempo trascorso. L’informazione può
anche essere letta nell’altro senso: se il problema è omogeneo e se si sa a priori
che:
u(x∗ , t∗ ) = max u(x, t) Con: (x∗ , t∗ ) ∈ QT
(x,t)∈Q̄T
Allora:
u(x, t) = cost ∀(x, t) ∈ QT
∂u
− µ∇2 u + v · ∇u + σu = f (x, t)
∂t
160
Ad esso, come sempre, vanno aggiunte la condizione iniziale di Cauchy e una delle
condizioni al bordo.
Per semplicità si consideri un problema monodimensionale omogeneo con con-
dizioni al contorno di Dirichlet omogenee:
u (x, t) − µuxx (x, t) + vux (x, t) + σu(x, t) = 0 (x, t) ∈ QT
t
u(0, t) = u(L, t) = 0 t ∈ (t0 , tf )
u(x, 0) = ϕ(x) x ∈ (0, L)
Uguagliando i coefficienti delle derivate con quelli presenti nel problema di diffu-
161
sione trasporto e reazione si ottiene il sistema:
v = −2µα α = − v
2µ
=⇒
σ = β − µα2 β = σ + v2
4µ
Per quanto riguarda le condizioni al contorno del problema in w(x, t) non ci sono
problemi, poiché:
w(0, t) = g(0, t)u(0, t) = 0
w(L, t) = g(L, t)u(L, t) = 0
∑
+∞ [ ( )2 ] ( )
v
x− µ nπ v2
+σ+ 4µ t nπ
u(x, t) = b̃n e 2µ L sin x
n=1
L
162
Capitolo 8
Problema di Laplace
Problema di Dirichlet
−∇2 u(x) = 0
x∈Ω
u(x) = g(x)
∂Ω
163
Problema di Neumann
−∇2 u(x) = 0
x∈Ω
∇u(x) · n̂ = g(x)
∂Ω
Problema di Robin
−∇2 u(x) = 0
x∈Ω
(∇u(x) · n̂ + γu(x)) = g(x)
∂Ω
164
Il primo membro dell’equazione può essere riscritto utilizzando il teorema della
divergenza:
∫ I ∫
( )
−µ ∇ · ∇u(x) dω = −µ ∇u(x) · n̂dσ = f (x)dω
Ω ∂Ω Ω
Da notare che questa condizione può essere interpretata come conservazione del
flusso di calore: tutto il calore generato internamente (funzione f (x)) coincide con
il calore smaltito sui bordi del dominio (funzione h(x)).
La funzione w(x) sarà soluzione del problema omogeneo (quindi di Laplace), poiché
le eventuali funzioni forzanti si semplificano nell’equazione differenziale. Inoltre,
165
lavorando sull’identità:
Condizioni di Dirichlet
u ∂Ω
=v ∂Ω
= g(x) =⇒ w ∂Ω
≡0
Condizioni di Robin
( ) ( )
−µ∇u· n̂−γu ∂Ω
= −µ∇v· n̂−γv ∂Ω
= g(x) =⇒ −µ∇w· n̂ ∂Ω
= γw
Condizioni di Neumann
∇u · n̂ ∂Ω
= ∇v · n̂ ∂Ω
= g(x) =⇒ ∇w · n̂ ∂Ω
≡0
166
E questo risultato implica che il secondo integrale è sempre nullo:
I
II = w∇w · n̂dσ = 0
∂Ω
∇w(x) ≡ 0 ∀x ∈ Ω
Essendo la derivata prima di w(x) sempre nulla, la funzione non può che essere
costante in tutto il dominio.
w = c = cost =⇒ u=v+c
u ∂Ω
=v ∂Ω
= g(x)
=⇒ c=0 =⇒ u≡v
167
che la derivata di w(x) è nulla, per cui si ritorna ad un problema di Dirichlet.
( )
− µ∇w · n̂ − γw ∂Ω
=w ∂Ω
=0
=⇒ c=0 =⇒ u≡v
∇w · n̂ ∂Ω
=0
=⇒ u = v + cost
168
Condizioni al contorno omogenee su tre lati
in cui le funzioni X(x) e Y (y) sono due funzioni a variabile reale del tipo:
X(x) : R −→ R Y (y) : R −→ R
Per prima cosa l’espressione della soluzione u(x, y) deve soddisfare l’equazione
differenziale, per cui, sostituendo si ottiene:
X ′′ Y ′′
X ′′ Y + XY ′′ = 0 =⇒ =− =k
X Y
in cui k non può essere che una costante reale, per via del fatto che entrambi i pro-
blemi sono definiti su una sola delle due variabili, che sono a loro volta indipendenti
tra loro.
Utilizzando il problema nella variabile x e le sue condizioni al contorno definite
al problema descritto dall’Eq. (8.1) si ottiene:
X ′′ (x) − kX(x) = 0
X(0) = X(L) = 0
Come si nota subito, la funzione g(x) è stata lasciata un attimo da parte, per ora
si considerano solamente condizioni al contorno omogenee.
169
Ancora una volta andando a risolvere il problema della funzione X(x) ci si deve
soffermare sul segno della costante k:
1 1
λL −λL
= e−λL − eλL ̸= 0 ∀L > 0
e e
c1 = c2 = 0
• costante nulla: k = 0
In questo caso la soluzione si riduce a un’espressione polinomiale, del tipo:
X(x) = c1 x + c2
170
Imponendo le condizioni al contorno si ottiene il seguente sistema:
X(0) = c = 0
2
X(L) = c L + c = 0
1 2
c1 = c2 = 0
In questo caso non è per nulla scontata l’unicità della soluzione, poiché il
determinante della matrice associata al sistema lineare è pari a:
1 0 nπ
= sin(λL) = 0 ⇐⇒ λ= n ∈ N∗
cos(λL) sin(λL) L
171
quello relativo alla variabile y:
Y ′′ (y) − ( nπ )2 Y (y) = 0
L
Y (H) = 0
c1,n e L H + c2,n e− L H = 0
nπ nπ
Yn (H) = 0 =⇒
nπ
c2,n = −c1,n e2 L H
[ nπ (y−H) nπ ]
nπ
H eL − e L (H−y)
= 2cn e L
2
[ − nπ (H−y) nπ ]
H −e + e L (H−y)
nπ L
= −2cn e L
2
nπ [ ]
= −2cn e L H sinh nπ L
(H − y)
[ nπ ]
= dn sinh L (H − y)
172
Ricomponendo tutti i pezzi si ottiene la seguente soluzione:
∑
+∞
[ nπ ]
u(x, y) = cn sin( nπ
L
x) sinh L
(H − y)
n=1
dove nei coefficienti cn sono state raggruppate tutte le varie costanti. A questo
punto occorre imporre che la soluzione soddisfi la condizione dettata dalla funzione
g(x):
∑
+∞
( nπ )
u(x, 0) = cn sin( nπ
L
x) sinh L
H = g(x)
n=1
Questo problema non è altro che la generalizzazione del caso precedente. In questa
nuova situazione le condizioni al contorno sono definite su ogni lato della regione
su cui è definita l’equazione; in tal caso il problema risulta essere:
∇2 u(x, y) = 0 (x, y) ∈ (0, L) × (0, H)
u = g i = 1, 2, 3, 4
Γi i
173
ratore differenziale del laplaciano (e quindi dell’equazione), utilizzando il principio
di sovrapposizione degli effetti; in tal modo si possono risolvere quattro problemi
più semplici, come quello analizzato nel caso precedente:
∇2 u(x, y) = 0 (x, y) ∈ Ω
per k = 1, 2, 3, 4
u = δ g i = 1, 2, 3, 4
Γi ik i
Importantissimo notare che ognuna delle quattro soluzioni è diversa dalle altre:
[ nπ ]
∑
+∞
sinh (H − y)
L[ ]
u1 (x, y) = b(1) nπ
n sin( L x) nπ
n=1
sinh L
H
[ nπ ]
∑
+∞
sinh (L − x)
H[ ]
u2 (x, y) = b(2) nπ
n sin( H y) nπ
n=1
sinh H
L
[ nπ ]
∑+∞
sinh y
u3 (x, y) = b(3) nπ
n sin( L x) [ nπL ]
n=1
sinh L H
[ nπ ]
∑
+∞
sinh x
u4 (x, y) = b(4) nπ
n sin( H y) [ nπ
H ]
n=1
sinh H L
174
solite cartesiane:
x = ρ sin(ϑ) ρ = x 2 + y 2
⇐⇒ ρ ∈ R+
0 , ϑ ∈ [−π, +π)
y = ρ sin(ϑ) ϑ = arctan y
x
∂2 1 ∂ 1 ∂2
∇2P = + +
∂ρ2 ρ ∂ρ ρ2 ∂ϑ2
A questo punto la soluzione sarà della forma (tecnica della separazione delle
variabili):
1 1
P′′ Θ + P′ Θ + 2 PΘ′′ = 0
ρ ρ
175
ricava:
P′′ (ρ) + ρ1 P′ (ρ) Θ′′ (ϑ)
− 1 = =k∈R
ρ2
P(ρ) Θ(ϑ)
Si parta dal problema più semplice, ossia quello in Θ(ϑ): al fine di di trovare una
soluzione accettabile per il problema, la funzione incognita deve per forza di cose
essere 2π-periodica (altrimenti dopo un giro risulta quantomeno discontinua).
• Costante positiva: k = λ2 > 0
Sotto questa ipotesi l’equazione differenziale diventa:
Quest’ultima è un’espressione che non può mai essere periodica, per cui la
soluzione è da scartare.
• Costante nulla: k = 0
Imponendo che la costante sia nulla l’equazione diventa:
Θ′′ (ϑ) = 0
Θ(ϑ) = c1 ϑ + c2
Θ(0) = Θ(2π) =⇒ c2 = c1 2π + c2
Θ0 (ϑ) = c0
176
• Costante negativa: k = −λ2 < 0
Da ultimo si impone che la constate sia negativa; da questo ne consegue che
l’equazione assume la forma:
Θ(−π) = Θ(π)
c1 cos(λπ) − c2 sin(λπ) = c1 cos(λπ) + c2 sin(λπ)
2c2 sin(λπ) = 0
Affinché non vi siano soluzioni banali, si deve per forza azzerare il seno; da
questo ne segue che:
λ = n ∈ N∗
Alla fine delle riflessioni relative alla funzione Θ(ϑ) si hanno due diverse tipo-
logie di funzioni accettabili:
1. costante nulla:
Il problema differenziale diventa:
1
P′′ (ρ) + P′ (ρ) = 0
ρ
177
Per risolvere questo nuovo problema non si fa altro che integrare:
P′′ (ρ) 1
′
=−
P (ρ) ρ
∫ ′′ ∫
P (ρ) 1
′
dρ = − dρ
P (ρ) ρ
ln|P′ (ρ)| = − ln|ρ| + k1
ρ
P(ρ) = c1 ln + c2
R
2. costante negativa:
Il problema differenziale a questo punto è completo:
178
sommatorie:
{( ) [ ]
ρ ∑
+∞
ρ
n
u(ρ, ϑ) = c0 + c̄0 ln + cn cos(nϑ) + dn sin(nϑ) +
R R
]}
n=1
( )−n [
ρ
+ c̄n cos(nϑ) + d¯n sin(nϑ)
R
Si consideri di essere nel caso della regione più semplice: una circolare. Applicando
le condizioni al contorno si ha che:
u(ρ, ϑ) ∂Ω
= g(ϑ)
∑
+∞
[ ]
g(ϑ) = a0 + an cos(nϑ) + bn sin(nϑ) ∀ϑ
n=1
A differenza delle altre volte l’integrazione per il calcolo dei coefficienti è legger-
mente diversa:
∫ +π ∫ ∫
1 1 +π 1 +π
a0 = g(ϑ)dϑ an = g(ϑ) cos(nϑ)dϑ bn = g(ϑ) sin(nϑ)dϑ
2π −π π −π π −π
179
Problema su una corona circolare
Questa volta non si sarà solamente un bordo, ma due (la circonferenza interna e
quella esterna):
{ }
Γ1 = (x, y) ∈ R 2 2
x +y = 2
R12
{ }
Γ2 = (x, y) ∈ R 2 2
x +y = 2
R22
Questa volta non si possono azzerare tutti i coefficienti delle funzioni che divergono
nell’origine, perché quest’ultima non appartiene al dominio. Andando a sviluppare
le funzioni g1 (ϑ) e g2 (ϑ) in serie di Fourier e imponendo le condizioni al contorno
180
si ottengono tre sistemi da risolvere per ricavare i coefficienti:
∫ +π
R1 1
c0 + c̄0 ln = g1 (ϑ)dϑ
R2 2π
∫ +π−π
1
c0 + c̄0 ln|1| = g2 (ϑ)dϑ
2π −π
( )n ( )−n ∫
R1 R1 1 +π
cn + c̄n = g1 (ϑ) cos(nϑ)dϑ
R2 R∫2 π −π
1 +π
cn 1n + c̄n 1−n = g2 (ϑ) cos(nϑ)dϑ
π −π
( )n ( )−n ∫
R 1 ¯ R1 1 +π
d n + dn = g1 (ϑ) sin(nϑ)dϑ
R2 R∫2 π −π
1 +π
dn 1n + d¯n 1−n = g2 (ϑ) sin(nϑ)dϑ
π −π
Tutti questi sistemi lineari risultano avere determinante non nullo se R1 ̸= R2 ,
quindi ammettono sempre una sola soluzione soluzione.
Formula di Poisson
181
Definizione 8.3.1 (Funzione armonica). Una funzione v(x) si dice armonica in
Ω ⊆ Rn se:
v(x) ∈ C 2 (Ω) ∇2 v(x) = 0
In altre parole, si può dire che una funzione è armonica se soddisfa l’equazione di
Laplace.
Dove |BR (x)| e |∂BR (x)| sono le misure dell’intorno: nella Tab. 8.1 sono riportati
i principali valori per diverse dimensioni.
182
Tabella 8.1: Principali valori delle misure di intorni sferici per diverse dimensioni.
Teorema 8.3.2 (Media integrale di funzioni armoniche). Data una funzione con-
tinua v(x) ∈ C 0 (Ω) che possiede le proprietà del Teo. 8.3.1 allora è anche una
funzione armonica e:
v(x) ∈ C ∞ (Ω)
Teorema 8.3.3 (Principio del massimo). Data una funzione v(x) armonica allora:
{ } { }
min v(x) ≤ v(x) ≤ max v(x)
x∈∂Ω x∈∂Ω
183
2. Si consideri il problema differenziale su Ω definito da:
∇2 u = 0
u =g
∂Ω
dove ug (x) è soluzione del problema associato. Siano dunque due funzioni
g1 e g2 due dati di Dirichlet al bordo e u1 e u2 le due soluzioni associate.
|u1 − u2 | ≤ max|g1 − g2 | ∀x ∈ Ω
x∈∂Ω
184
Parte III
185
Capitolo 9
In particolare si ha che: ∫
gn (x)dx = 1 ∀n
R
Il delta di Dirac è quindi la funzione impulso, cioè vale zero sempre tranne nell’o-
rigine in cui assume il valore uno. In realtà il punto su cui avviene l’impulso può
186
essere spostato in un punto c qualsiasi, tramite le seguenti notazioni.
δc (x) δ(x − c)
La derivata della funzione di Heaviside può essere identificata nel delta di Dirac:
d
H(x − c) = δc (x)
dx
187
Si nota immediatamente che la funzione f (x) = δ 1 (x) non è una funzione continua
2
nel dominio, quindi non può esistere la formulazione forte al problema appena
presentato; tuttavia nella pratica comune si può avere esperienza tutti i giorni di
questo fenomeno.
Non esiste quindi una formulazione forte al problema, ma può esistere una
formulazione debole. Per prima cosa si può integrare da entrambi i lati:
∫ ∫ c
x x x ∈ (0, 1/2)
−u′ (x) = − u′′ (τ )dτ = δ 1 (τ )dτ + c =
0 0
2 c + 1 x ∈ [1/2, 1)
−c x ∈ (0, 1/2)
′
u (x) =
−(c + 1) x ∈ [1/2, 1)
che corrisponde alla soluzione del problema. Quest’ultima è una funzione continua,
ma non derivabile neanche al primo ordine; per questo motivo non si può dire che
è una soluzione al problema posto in formulazione forte.
La tecnica appena vista prende il nome di formulazione del problema diffe-
renziale debole, ossia si cerca di integrare una funzione, eventualmente sfruttando
alcune particolari proprietà di integrazione
188
Capitolo 10
Integrale di Lebesgue
Cioè si utilizza una funzione g di poco più regolare di f (e quindi più facile da
integrare) per svolgere l’operazione di integrazione.
Quindi si ha che: la funzione f è integrabile se e solo se la funzione g risulta
integrabile. Inoltre si hanno le seguenti proprietà:
189
Derivazione nel senso delle distribuzioni
∀ v, u ∈ V v+u=w ∈V
∀v ∈ V α∈R αv = w ∈ V
1. Elemento neutro:
v+0=v 1v = v
2. Elemento opposto:
1
v + (−v) = 0 v =1
v
3. Commutativa:
v+u=u+v αv = vα
u + (v + w) = (u + v) + w (αβ)v = α(βv)
190
5. Distributiva.
1. l’insieme R;
2. i vettori di Rn ;
3. i polinomi di grado n, Pn ;
5. l’insieme delle funzioni continue su [a, b], intervallo generico, C 0 ([a, b]).
1. Positività:
2. Omogeneità:
∀λ ∈ R ∀v ∈ V ∥λv∥ = |λ|∥v∥
3. Disuguaglianza triangolare
∀ v, u ∈ V ∥x + y∥ ≤ ∥v∥ + ∥u∥
d(x, y) = ∥x − y∥
191
In questo caso si dice che lo spazio vettoriale V è normato, cioè dotato di norma,
e metrico, poiché su di esso è definita una distanza.
d(xn , xm ) = ∥xn − xm ∥ −→ 0+ n, m −→ +∞
Non è assolutamente garantito il contrario, cioè non è detto che una sequenza
di Cauchy converga a x.
∥·∥ = |·|
192
conclude dicendo che la successione scelta è una sequenza di Cauchy, ma non è
una sequenza di elementi in X = Q
su tutti i punti x ∈ Ω, fatta eccezione di alcuni punti isolati in cui la funzione può
essere singolare.
Sulla questione della norma bisogna in genere stare molto attenti: ad esempio,
si consideri lo spazio delle funzioni continue C 0 (Ω). Tale spazio è dotato di una
norma (quindi è normato), ma non è completo se si usa come norma quella di
L1 (Ω) o quella di L2 (Ω), ad esempio la funzione f (x) = 1 presenta una forma in-
determinata quando p = +∞. Tuttavia lo spazio C 0 (Ω) diventa normato completo
quando si scegli un’altra tipologia di norma: quella del massimo:
193
Teorema 10.1.3 (Disuguaglianza di Hölder). La disuguaglianza di Hölder riguar-
da la norma del prodotto di due funzioni:
∫
1 1
|f g|dω ≤ ∥f ∥Lp ∥g∥Lq con: + = 1 [1 ≤ p, q < +∞)
Ω p q
⟨·, ·⟩ : X × X −→ R
2. commutativa:
∀ x, y ∈ X ⟨x, y⟩ = ⟨y, x⟩
3. distributiva:
194
Gli spazi di Hilbert son spazi funzionali necessari per i metodi matematici
avanzati. Da notare che non tutti gli spazi vettoriali di Banac (Def. 10.1.7) Lp (Ω)
non sono spazi di Hilbert.
Ossia si ha una sola proiezione ortogonale del vettore x ∈ H su V data dal punto
Pv x tale da minimizzare la distanza tra x e il generico vettore v ∈ V .
x∈V ⇐⇒ Pv x = x
195
Cioè lo spazio funzionale in cui tutte le funzioni in esso contenute e le loro derivate
sono elementi di L2 (Ω). Più in generale si dice che:
{ }
k
H (Ω) = v ∈ L (Ω)
2
D v ∈ L (Ω) ∀ α
α 2
|α| ≤ k
n
H k (Ω) ⊂ Cm (Ω̄) se: k ≥ m +
2
Quindi si sa per certo che tutte le funzioni in una dimensione contenute nello spazio
di Sobolev H 1 (Ω) sono continue. Questo non vale se le dimensioni aumentano
n ≥ 2 e lo spazio non cambia:
n
H k (Ω) ⊂ C 0 (Ω) se: k > ≥1
2
196
Teorema 10.1.7 (Disuguaglianza di traccia). Assegnato v ∈ H 1 (Ω) dove il bordo
∂Ω risulta sufficientemente regolare allora esiste un’applicazione lineare:
γ0 H 1 (Ω) −→ L2 (Ω)
Tale che:
γ0 v = v ∂Ω
e ∥γ0 v∥L2 ≤ c∥v∥H 1
Definizione 10.1.13 (Norme equivalenti). Due norme ∥·∥ e |·| si dicono equiva-
lenti se esistono due costanti positive c1 , c2 ∈ R+ tali che:
197
Da cui segue che:
|v|2H 1 ≤ ∥v∥2L2 + |v|2H 1 = ∥v∥2H 1
∥v∥L2 ≤ CΩ |v|H 1
Bene adesso che si è maggiorata la funzione si calcola la sua norma su L2 (a, b):
∫ b ∫ b (∫ x )
′ 2
∥v∥2L2 = 2
v (x)dx = (x − a) |v (x)| ds dx
a a a
b−a
CΩ = √
2
198
10.2 Funzionali e forme
I funzioni e le forme sono degli elementi essenziali per riscrivere i problemi diffe-
renziali in forma debole. Di fatto essi rappresentano una sorta di generalizzazione
delle classiche funzioni.
F : H −→ R
I funzionali possono avere o non avere alcune proprietà, a seconda di come sono
definiti. In particolare, possono essere:
• Lineari:
F (u + λv) = F (u) + λF (v) ∀ u, v ∈ H, λ ∈ R
• Limitati:
∃c > 0 |F (v)| < c∥v∥H ∀v ∈ H
• Continui:
F : H −→ R
|F (v)|
∥F ∥H ∗ = sup
v∈H ∥v∥H
v̸=0
199
Teorema 10.2.2 (Teorema di rappresentazione di Riesz). Sia H uno spazio di
Hilbert, allora ∀ F ∈ H ∗ :
a : H × H −→ R
Anche le forme, come i funzionali, possono avere o non avere alcune importanti
proprietà; tra esse si riportano:
• Bilinearità:
a(u, λv + µz) = λa(u, v) + µa(u, z)
∀ u, v, z ∈ H, ∀ λ, µ ∈ R
a(λu + µz, v) = λa(u, v) + µa(z, v)
• Simmetria:
a(u, v) = a(v, u) ∀ u, v ∈ H
• Continuità:
• Coercività:
∃α > 0 a(u, u) ≥ α∥u∥H ∀u ∈ H
200
Capitolo 11
Come già accennato all’interno della Sez. 9 a pagina 186 molti dei problemi dif-
ferenziali più complessi non ammettono una soluzione classica, ma potrebbero
ammetterne una in forma debole.
a(u, v) = F (v) ∀v ∈ V
allora esiste una sola soluzione al problema posto in forma debole u e inoltre:
1
∥u∥V ≤ ∥F ∥V ∗
α
201
Dimostrazione 11.0.1 (Norma della soluzione). Andando a maggiorare il fun-
zionale, tramite la definizione della norma dello spazio duale (Def. 10.2.2), e
minorando la forma bilineare tramite la proprietà di coercività (Def. 10.2.3), si
ottiene:
α∥u∥2V ≤ a(u, u) = F (u) ≤ ∥F ∥V ∗ ∥u∥V
α∥u∥2V ≤ ∥F ∥V ∗ ∥u∥V
1
∥u∥V ≤ ∥F ∥V ∗
α
u J(v) v=u
= min J(v)
v∈V
Tra tutte le funzioni possibili si sceglie la funzione v = δ. In tal modo il valore del
202
funzionale in ṽ = u + δ è pari a:
∫ L
1
J(u + δ) = J(u) + µ0 (δ ′ )2 dx ≥ J(u)
2 0
Questo perché il termine della parentesi quadra è pari a zero, essendo soluzione.
Dato che si sa per certo che la funzione u è il minimo dell’energia si può scrivere:
∂J
(u + δw) = 0 ∀w ∈ V
∂δ
Quindi per forza di cose la funzione u è soluzione del problema variazionale astrat-
to.
203
11.1.1 Problema di Poisson-Dirichlet
Condizioni al contorno omogenee
Riprendendo l’Eq. (11.1) e integrando per parti il temine a sinistra per parti si
ottiene: ∫ L ∫ L
L
′ ′ ′
−(µ0 u v) + µ0 u v dx = f vdx ∀v ∈ V
0 0 0
Avendo scelto V = H01 (0, L) il primo termine è sempre nullo, poiché la ogni funzione
v si annulla in corrispondenza degli estremi del dominio. In tal modo si ottiene
dunque: ∫ ∫
L L
µ0 u′ v ′ dx = f vdx
0 0
Per la disuguaglianza di Hölder (Teo. 10.1.3), entrambi gli integrali sono definiti
se:
u, u′ , v, v ′ ∈ L2 (0, L)
204
cosa che va benissimo, poiché si è appena scelto lo spazio H01 (0, L) per queste
funzioni. Il problema variazionale astratto (o problema in forma debole) risulta
così definito:
∫ L ∫ L
′ ′
u∈V µ0 u v dx = f vdx ∀v ∈ V con V = H01 (0, L)
|0 {z } | {z }
0
a(u,v) F (v)
Si evince immediatamente che è lineare sia in u che in v,e quindi bilineare, poiché
le operazioni di derivazione e integrazione sono entrambe lineari. La forma è anche
continua, poiché per la disuguaglianza di Schwarz:
∫ L
|a(u, v)| ≤ |µo u′ v ′ |dx ≤ µ0 ∥u′ ∥L2 ∥v ′ ∥L2
0
Questa non è ancora la definizione di forma continua per che le norme sono en-
trambe calcolate sullo spazio L2 , che non è lo spazio V scelto. Tuttavia, proprio
per come stata definita la norma su H 1 si ha che:
√
∥u∥H 1 = ∥u∥2L2 + ∥u′ ∥2L2 ≥ ∥u′ ∥L2
Lo stesso passaggio può essere riproposto per la funzione v. Risulta quindi che:
205
Dato che lo spazio scelto è V = H01 allora è possibile utilizzare la disuguaglianza
di Poincaré (Teo. 10.1.9). Per questo motivo è lecito scrivere:
206
bordo Rg , tale per cui:
x
Rg (x) = g1 + (g2 − g1 )
L
A questo punto si può riscrivere la funzione soluzione come somma della funzione
Rg e della funzione ū, ancora incognita:
Ancora una volta il secondo termine dell’integrazione per parti si annulla, essendo
la funzione test v sempre nulla al bordo. Riscrivendo i termini rimanenti si ha che:
∫ L ∫ L
′ ′
µ0 ū v dx = (f v − µ0 Rg′ v ′ )dx
| 0
{z } |0 {z }
a(ū,v) F (v)
Come si evince facilmente la forma bilineare non è cambiata, quindi sarà sempre
continua e coerciva. Il problema è capire se il funzionale F (v), in cui sono rientrate
207
anche le condizioni al contorno risulta ancora lineare e continuo.
Sulla prima richiesta non ci sono problemi, le operazioni di derivazione, som-
ma algebrica e integrazione sono lineari, quindi la combinazione di queste tre è
anch’essa lineare. Rimane dunque dimostrare la continuità:
∫ L ∫ L
|F (v)| ≤ |f v|dx + µ0 Rg′ v ′ dx
0 0
= ∥f v∥L1 + µ0 Rg′ v ′ L1
≤ ∥f ∥L2 ∥v∥L2 + µ0 Rg′ L2 ∥v ′ ∥L2
≤ ∥f ∥L2 ∥v∥H 1 + µ0 Rg′ L2 ∥v∥H 1
( )
= ∥f ∥L2 + µ0 Rg′ L2 ∥v∥H 1
c = ∥f ∥L2 + µ0 Rg′ L2
Quindi per il teorema di Lax-Milgram (Teo. 11.0.2) esiste una sola soluzione ū ∈ V
al problema variazionale astratto. Naturalmente una volta risolto il problema in ū,
si deve sommare la soluzione di rilevamento Rg per trovare la soluzione al problema
principale.
208
Integrando il primo termine per parti si ha che:
L ∫ L ∫ L
′ ′ ′
−µ0 u (x)v(x) + µ0 u v dx = f vdx
0 0 0
Questa volta i termini valutati non spariscono, poiché non è stata fatta ancora nes-
suna ipotesi sullo spazio funzionale del problema V , ne sulla forma della funzione
test v. Non c’è altra scelta che valutare i termini in corrispondenza del bordo:
∫ L ∫ L
′ ′ ′ ′
− µ0 u (L) v(L) + µ0 u (0) v(0) + µ0 u v dx = f vdx
| {z } | {z } 0 0
h2 h2
Non a caso le condizioni al contorno sono finite nel funzionale F (v). Questo succede
sempre se le condizioni al contorno non sono quelle di Dirichlet omogenee.
Dato che le condizioni al bordo di Neumann, Robin, miste e Dirichlet non
omogenee sono talmente importanti che finiscono all’interno della formulazione
debole, sono anche dette condizioni naturali, mentre quelle di Dirichlet omogenee
sono chiamate condizioni essenziali, perché la loro presenza ha effetto solo sulla
scelta dello spazio funzionale adeguato, ma non rientrano nel problema principale.
Manca comunque da definire un spazio funzionale in cui ambientare il problema
e da verificare che tutte le ipotesi di Lax-Milgram siano ancora verificate. Come
si nota dalla condizioni al contorno, non per forza la soluzione u deve annullarsi
in corrispondenza degli estremi e quindi non è più possibile scegliere lo spazio H01 .
Nonostante questa complicazione, rimane più che lecito scegliere lo spazio H 1 . Con
questo tipo di scelta il problema variazionale diventa:
209
Si passa ora alla verifica delle ipotesi di Lax-Milgram. La forma è sempre la stessa,
quindi non ci sono problemi. Il funzionale invece è cambiato, quindi va rivisto. La
linearità è al solito garantita, mentre la continuità va verificata:
∫ L
|F (v)| ≤ |f v|dx + |h2 ||v(l)| + |h1 ||v(0)|
0
Anche in questo caso tutte le ipotesi sono rispettate, quindi esiste solo una solu-
zione al problema in forma debole.
210
Capitolo 12
Metodo di Galërkin e
approssimazione ad elementi finiti
211
base:
∑
+∞
{ }+∞
w∈V ⇐⇒ w= wj ηj con: wj i=1
∈R
j=1
Vh ⊂ V dim(Vh ) = Hh < +∞
A questo punto si può costruire una sequenza di spazi, ognuno avente dimensione
superiore a tutti quelli precedenti e dimensione inferiore a quelli successivi:
dove ovviamente h1 < h2 < · · · < hNh . Inoltre, proprio per come è definita la
sequenza si ha che:
lim Nh = +∞ lim Vh = V
h→0 h→0
Tutti i sottospazi sono sottospazi vettoriali, quindi ogni elemento ad essi appar-
tenente può ancora essere scritto come combinazione lineare degli elementi della
base:
∑Nh
vh ∈ Vh ⇐⇒ vh = vj φj
j=1
in cui le φj sono le funzioni di base, mentre i pesi vj sono coefficienti reali detti
variabili di controllo.
a(uh , vh ) = F (vh ) ∀ vh ∈ Vh
212
colare, in funzione della scelta della base {φj }j=1
Nh
di Vh si avrà lo specifico metodo
di Galërkin, in cui le incognite saranno {uj }j=1 . Il problema variazionale astratto
Nh
a(uh , vh ) = F (vh ) ∀ vh ∈ Vh
a(uh , φi ) = F (φi ) i = 1, 2, . . . , Nh
(∑
Nh )
a uj φj , φi = F (φi ) i = 1, 2, . . . , Nh
j=1
A questo punto si sfrutta la proprietà di bilinearità della forma a(·, ·); in tal caso
il problema può essere riscritto come:
∑
Nh
a(φj , φi ) uj = F (φi ) i = 1, 2, . . . , Nh
j=1
| {z } | {z }
[A]ij fi
Au = f
213
In poche parole, si è appena dato prova che uh ∈ Vh ⊂ V è soluzione del problema
variazionale astratto se e solo se e soluzione del problema di Galërkin.
Le proprietà della matrice A, chiamata matrice di rigidezza o matrice di stiff-
ness, dipendono dalla forma a(·, ·) associata e dalla base che viene adottata per
risolvere il problema.
∑
Nh
[ ] ∑
Nh
T
v Av = vi A ij vj = vi a(φj , φi )vj
i,j=1 i,j=1
214
Esistenza e unicità
Stabilità
1
∥uh ∥Vh ≤ ∥F ∥V ∗
α h
Da notare che è la stessa dimostrazione proposta per l’ultimo punto del teorema
di Lax-Milgram.
215
Convergenza
a(u − uh , vh ) = 0 ∀ vh ∈ Vh
Dato che la prima equazione è più generale, per via del fatto che Vh ⊂ V , si può
scegliere una particolare funzione test v = vh ∈ Vh . Il tutto può essere riscritto nel
sistema:
a(u, v ) = F (v )
h h
∀ vh ∈ Vh
a(u , v ) = F (v )
h h h
216
Inoltre, è bene sottolineare che se la forma a(·, ·) è simmetrica, allora è possibile
introdurre la norma dell’energia come:
√
∥u∥a = a(u, u)
Quindi nel caso in cui il dominio della forma sia V × V = R2 , si può osservare
che il Teo. 12.1.4, scritto attraverso la norma dell’energia, implichi che il prodotto
scalare fra la differenza u − uh e una qualsiasi funzione test vh ∈ Vh sia sempre
nullo.
Infatti, la soluzione al metodi di Galërkin uh non è altro che la proiezione della
soluzione u sul sottopsazio Vh ⊂ V (ortogonale nel senso della prodotto scalare
indotto dalla forma); quindi la soluzione al problema di Galërkin è certamente
la migliore approssimazione della soluzione u del problema in forma debole, una
volta scelto lo spazio Vh . Questo concetto può essere riassunto dalla proprietà della
consistenza forte:
u ∈ Vh =⇒ uh ≡ u
M
∥u − uh ∥V ≤ inf ∥u − vh ∥V
α vh ∈Vh
a(u − uh , u − uh ) = a(u − uh , u − vh + vh − uh )
217
forma continua e coerciva, si ha che:
α∥u − uh ∥V ≤ M ∥u − vh ∥V ∀ vh ∈ Vh
M
∥u − uh ∥V ≤ inf ∥u − vh ∥V
α vh ∈Vh
Per h → 0 =⇒ Vh → V
( )
=⇒ lim inf ∥v − wh ∥V = 0 ∀v ∈ V
h→0 wh ∈Vh
218
scelta particolare dello spazio Vh .
Per prima cosa si partiziona l’insieme Ω̄ in N + 1 sottointervalli, chiamati
elementi di griglia:
Kj = (xj−1 , xj ) j = 1, 2, . . . , N + 1
a = x0 , x1 , x2 , . . . , xN , xN +1 = b
Non è detto che tutti questi nodi siano equispaziati. Nel caso non lo fossero il
parametro h è pari al massimo di tutte le distanze fra nodi consecutivi:
h= max hj hj = xj − xj−1
j=1,...,N +1
Definizione 12.2.1 (Spazio elementi finiti). Si chiama spazio degli elementi finiti
di grado r ≥ 1 lo spazio funzionale:
{ }
Xhr (Ω) = vh ∈ C 0 (Ω̄) vh Kj
∈ Pr ∀Kj ∈ τh
Dato che in generale è utile non complicarsi troppo la vita, si usano essenzial-
mente due diverse tipologie di elementi finiti: lineari e quadratici.
Il passo successivo è quello di costruire le basi {φj }N
j=0 dello spazio a elementi
h
219
Si avrà quindi che la base {φj }N +1
j=0 sarà composta da N +2 funzioni, definite come:
x −x
1 x ∈ [x0 , x1 ]
φ0 = x1 − x0
0 x ̸∈ [x0 , x1 ]
x − xi−1
x ∈ [xi−1 , xi )
xi − xi−1
xi+1 − x
φi = x ∈ [xi , xi+1 ] i = 1, 2, . . . , N
xi+1 − xi
0 x ̸∈ [xi−1 , xi+1 ]
x−x
N
x ∈ [xN , xN +1 ]
φN +1 = xN +1 − xN
0 x ̸∈ [xN , xN +1 ]
In generale, si tende a ricercare delle basi a elementi finiti con supporto (cioè la
distanza massima tale che φj φi ̸= 0) più piccolo possibile, in modo che la matrice
di rigidezza riempia di zeri e quindi sia più facile da risolvere.
La base che è appena stata ottenuta è detta base lagrangiana (Def. 3.1.2 a
pagina 50), ed è la stessa che si è trovata nel capitolo sull’approssimazione di
funzioni. Come già analizzato in precedenza, questa tipologia di polinomi ha la
proprietà di assumere valore unitario quando vengono valutati nel nodi in cui sono
centrati e di esser nulli in tutti gli altri nodi.
φj (xi ) = δij
dove δij è il delta di Kronecker (Def. 3.1.3). Inoltre tutte le incognite uh possono
essere riscritte come come combinazione lineare degli elementi della base Inoltre,
220
φ0 φi
1 1
x x
x0 x1 xN +1 xi−1 xi xi+1
(a) Primo elemento di una base generica (b) Elemento i-esimo di una base generi-
a elementi finiti lineari monodimensiona- ca a elementi finiti lineari monodimensio-
le. nale.
φN +1 φ
1 1
x x
x0 xN xN +1 x0 xN +1
(c) Ultimo elemento di una base generica (d) Intera base a elementi finiti lineari
a elementi finiti lineari monodimensiona- monodimensionale.
le.
221
proprio per la proprietà del delta di Kronecker si ha che:
∑
N +1 ∑
N +1
uh (x) = uj φj (x) =⇒ uh (xi ) = uj φj (xi ) = ui
j=0 j=0
x = ϕi (ξ) = xi + ξ(xi+1 − xi )
In cui ϕi (x) : [0, 1] −→ [xi , xi+1 ] è detta funzione di traccia e serve, una volta
invertia, per passare da una mappa sempre uguale alla generica funzione base
φi (x). Definendo due funzioni φ̂0 (ξ) e φ̂1 (ξ) con i valori di ξ ∈ [0, 1], si ottiene la
mappa:
φ̂ (ξ) = 1 − ξ
0
φ̂ (ξ) = ξ
1
A questo punto si può invertire la funzione di traccia, per passare alle funzioni di
base, cioè ∀ i = 0, 1, . . . , N :
( )
φi (x) = φ̂0 ξ(x)
( )
φi+1 (x) = φ̂1 ξ(x)
x − xi
con: ξ(x) = ϕ−1i (x) =
xi+1 − xi
Notare le due funzioni corrispondono a metà della funzione di base centrata sul
nodo sinistro e metà della funzione di base che è centrata sul nodo destro; nella
Fig. 12.2 è riportata la mappatura di un elemento generico della base.
222
φ̂0 , φ̂1 φi , φi+1
1 1
ξ x
1 xi xi+1
Figura 12.2: Mappatura per una base a elementi finiti lineari in una dimensione.
Quindi in questo caso non si hanno più delle rette, ma delle parabole. Ovviamen-
te data la maggiore complessità del polinomio sono necessari tre punti per ogni
parabola.
In altre parole si fa in modo che i nodi di indice pari siano i vertici delle parabole
rivolte verso il basso, mentre quelli di indice dispari degli zeri di altre parabole
223
φpari φdisp
1 1
x x
xi−2 xi−1 xi xi+1 xi+2 xi−1 xi xi+1
(a) Elemento finito quadratico pari. (b) Elemento finito quadratico dispari.
x
x0 x2N +2
φj (xi ) = δij
Questa tipologia elementi finiti è riportata nella Fig. 12.3; in particolare si nota
che il supporto degli elementi quadratici pari è il doppio di quelli dispari.
Anche in questo caso è possibile adottare la tecnica della mappatura della base,
224
che risulta un po’ più complicata di quella lineare. Nello specifico:
φ̂ = (1 − ξ)(1 − 2ξ)
0 ( )
φ̂1 = 4ξ(1 − ξ) =⇒ φi (x) = φ̂k ξ(x)
φ̂ = ξ(2ξ − 1)
2
uh ∈ Vh ⊂ V a(uh , vh ) = F (vh ) ∀ vh ∈ Vh
225
φ
uh (0) uh (L)
x
L
Figura 12.4: Base degli elementi finiti lineari per il problema di Poisson con condizioni
di Dirichlet omogenee; in azzurro sono segnati tutti gli elementi della base presenti, in
viola quelli automaticamente esclusi dalle condizioni al contorno.
Si analizzi per prima il caso più semplice, quello a elementi finiti lineari. Si impone
quindi che:
{ }
Vh = Xh1 = vh ∈ C (Ω̄)
0
vh Kj
∈ P1 ∀ Kj ∈ τh
Notare che la scelta dello spazio (già a livello del problema in forma debole) era
stata V = H01 , quindi lo spazio adesso risulta essere:
In tal modo si sa per certo che la funzione deve essere nulla agli estremi, per cui si
può evitare già dal principio di costruire la base ad elementi finiti in quelle zone;
in Fig. 12.4 è riportata la base scelta.
Adesso è necessario costruire la soluzione, con la solita combinazione lineare
226
(da notare gli estremi della sommatoria!):
∑
N
uh = uj φj (x)
j=1
A questo punto si hanno tutti gli strumenti per costruire la matrice di rigidezza
A e il vettore delle forzanti f . Per prima cosa è bene fare una riflessione sulle loro
dimensioni:
dim(Vh ) = N =⇒ A ∈ RN ×N , f ∈ RN
∑
N N (∫
∑ L ) ∫ L
a(φj , φi )uj = µφ′j φ′i dxuj = f φi dx i = 1, 2, . . . , N
j=1 j=1 0 0
Per prima cosa, si assembla la matrice di rigidezza. Per farlo, il primo passo è
calcolare le derivate degli elementi finiti.
Essendo loro delle rette, le derivate prime non saranno alto che i coefficienti
angolari (cfr. Fig 12.5):
1 x ∈ (xi−1 , xi )
φ′i (x) = h
− 1 x ∈ (xi , xi+1 )
h
• Se |i − j| > 1 allora per forza almeno una delle due derivate (j o i che sia)
si annulla, quindi:
227
φi (x) φ′i (x)
1
1 h
x
xi−1 xi xi+1
x
− h1
xi−1 xi xi+1
quindi: ∫ ( )
L
1 1 1
Aij = µ − dx = −µ
0 h h h
228
gli elementi poi saranno moltiplicati per un fattore µ/h. In formule:
2 −1 0 0 ...
..
−1 2 −1 .
..
µ 0 −1 2 .
A=
h .. ..
. .
..
. 2 −1
0 ... −1 2
Quanto riguarda il vettore delle forzanti il discorso è più semplice: ogni com-
ponente è il risultato di un’integrazione definita:
∫ L
fi = f (x)φi (x)dx i = 1, 2, . . . , N
0
In questo caso si vogliono utilizzare gli elementi finiti quadratici, al posto di quelli
lineari. Per questo motivo si ha che:
229
φ
uh (0) uh (L)
x
L
Figura 12.6: Base degli elementi finiti quadratici per il problema di Poisson con condi-
zioni di Dirichlet omogenee; in azzurro sono segnati tutti gli elementi della base presenti,
in viola quelli automaticamente esclusi dalle condizioni al contorno.
Questo perché gli elementi finiti quadratici hanno un supporto pari a 4h quando
sono dispari, mentre se sono pari hanno un supporto pari a 2h.
In ogni caso la matrice di rigidezza A ha una serie di importanti proprietà che
val la pena di sottolineare:
1. A è una matrice definita per bande, cioè composta da elementi non nulli a
bande diagonali;
c
K2 (A) = c>0
h2
230
zione, poiché Ks (A) → +∞ (tutto questo senza neanche contare il costo
computazionale. . . ).
Il concetto alla base è pressoché identico. Per prima cosa si crea la base di elementi
finiti; dato un insieme Ω definito su più dimensioni, si ha che:
{ }
Vhr (Ω) = vh ∈ C 0 (Ω̄) vh Kj
∈ Pr ∀ Kj ∈ τh
231
Dove lo spazio funzionale scelto Hd1 (Ω) prevede che tutte le funzioni v abbiano uno
zero in corrispondenza del bordo sinistro:
{ }
Hd1 (Ω) = v ∈ H 1 (Ω) v(0) = 0
Dato che nell’estremo sinistro la soluzione presenta uno zero, è possibile tralasciare
il primo elemento della base ad elementi finiti, la quale quindi è esprimibile con la
notazione {φi (x)}N +1
i=1 (cfr. Fig 12.7).
Da ultimo, si procede con l’assemblaggio della matrice A e del vettore delle
forzanti f . In realtà, questi sono molto simili al problema con le condizioni al
contorno di Dirichlet omogenee (Par. 12.3.1), fatta eccezione per le dimensioni e
gli ultimi elementi:
A ∈ R(N +1)×(N +1) f ∈ RN +1
La matrice di rigidezza A è del tutto analoga a quella precedente, tranne che per
232
φ
uh (0)
x
L
Figura 12.7: Base degli elementi finiti per il problema di Poisson misto con condizioni
di Dirichlet omogenee a sinistra e di Robin a destra; in azzurro sono segnati tutti gli
elementi della base presenti, in viola quello automaticamente esclusi dalle condizioni
al contorno di Dirichlet. Sulla destra è incluso quello in verde dove sono applicate le
condizioni di Robin, che non per forza implicano che la funzione uh presenti uno zero in
L.
233
Per quanto riguarda il vettore delle forzanti f il discorso è molto più facile:
∫L
f (x)φ0 (x)dx
0
∫ 0
L
∫ 0
f (x)φ 1 (x)dx
0
L f (x)φ2 (x)dx 0
0
f= .. + .
.
..
∫L
f (x)φ N (x)dx
0
∫
L
0
0
f (x)φ N +1 (x)dx q 2
Infatti il vettore è del tutto analogo a quello precedente, fatta eccezione per l’ultimo
termine, che vede invece l’aggiunta delle condizioni al contorno di Robin.
µ, σ ∈ R+ β∈R
234
Dove la forma e il funzionale valgono rispettivamente:
∫ L ∫ L ∫ L
′ ′ ′
a(u, v) = µ u v dx + β u vdx + σ uvdx
0 0 0
∫ L
F (v) = f vdx + q1 v(0) + q2 v(L)
0
Da notare che lo spazio vettoriale non ha più gli zeri agli estremi. Per questioni di
semplicità si opta per uno spazio di elementi lineare, ottenuto da una partizione
di N + 2 nodi equispaziati, con passo h:
{ }
L
Vh = Xh1 = vh ∈ C 0 (Ω̄) vh Kj
∈ P1 ∀ Kj ∈ τ h h=
N +1
Al contrario del problema di Dirichlet, questa volta bisogna per forza di cose
utilizzare una base ad elementi finiti completa (Fig. 12.8). Di conseguenza, le
dimensioni della matrice e dei vettori sono:
∑
N +1
A ∈ R(N +2)×(N +2) f ∈ RN +2 uh (x) = uj φj (x)
j=0
Per prima cosa si parte dalla forzante f , più semplice da calcolare. Dal proble-
ma, si ha che:
∫ L
fi = f (x)φi (x)dx + q1 φi (0) + q2 φi (L) i = 0, 1, . . . , N + 1
0
Dato che gli ultimi due termini danno un contributo solo sul primo e sull’ultimo
235
φ
x
L
Figura 12.8: Base degli elementi finiti per il problema di diffusione reazione e trasporto
stazionario con condizioni di Neumann; sono evidenziati in viola tutti gli elementi degli
estremi, presenti anch’essi nel problema.
A=K+B+M
236
elementi [K]11 e [K]Nh Nh :
∫ L ( )2 µ
[K]11 = µ φ′0 dx = = [K]Nh Nh
0 h
In questo caso non si ha più il coefficiente due perché si sta integrando sul primo
e l’ultimo elemento finito, che sono pari alla metà di uno standard. La matrice K
risulta quindi essere:
1 −1 0 0 ...
..
−1 2 −1 .
..
µ 0 −1 2 .
K= det(K) = 0
h .. ..
. .
..
. 2 −1
0 ... −1 1
Per assemblare la matrice B, bisogna distinguere quattro casi diversi (cfr. Fig. 12.9):
• Se |i − j| > 1 allora per forza almeno una fra la derivata o la funzione base
si annulla, quindi:
237
φ
1
h
≈
xi−1 xi xi xi+1 xi−1 xi xi+1
− h1
Figura 12.9: Base degli elementi finiti per il problema di diffusione reazione e trasporto
stazionario per il calcolo della matrice B. In ordine sono rappresentati il caso j = i − 1,
j = i + 1 e infine j = i.
zata in due parti, nella prima la derivata risulta positiva, mentre la seconda
sarà negativa. I due integrali avranno quindi lo stesso valore, ma con segno
opposto:
∫ L ∫ xi+1 ∫ xi ( )
1 1 1h 1h
Bij = β φ′j φi dx =β φi dx − β φi dx = β − =0
0 xi h xi−1 h h2 h2
• Se |i − j| > 1 allora almeno una delle due funzioni risulta esser nulla, quindi:
238
• Se j − i = ±1 allora una funzione è crescete, mentre l’altra è decrescente.
Ne risulta quindi che:
∫ L
σh
Mij = σ φj (x)φi (x)dx =
0 6
Anche in questo caso si può passare alla formulazione debole del problema, la quale
è sintetizzabile come:
239
In cui la forma e il funzionale sono definiti esattamente come nel caso precedente,
fatta eccezione per le funzioni µ(x), β(x) e σ(x), che questa volta non possono
uscire dall’integrale:
∫ L ∫ L ∫ L
′ ′ ′
a(u, v) = µu v dx + βu vdx + σuvdx
0 0 0
∫ L
F (v) = f vdx
0
A questo punto si hanno tutti gli strumenti per passare al metodo di Galërkin-
elementi finiti; la matrice di rigidezza risulta quindi:
∫ ∫ ∫
[ ] L L L
A ij = µ(x)φ′j φ′i dx + β(x)φ′j φi dx + σ(x)φj φi dx
0 0 0
∑
N +1
fi ≃ ᾱq f (ȳq )φi (ȳq )
q=0
240
12.5 Elementi finiti per un problema di diffusio-
ne
Il caso del problema di diffusione è notevolmente più complicato di tutti gli altri
riscontrati finora: infatti questa è la prima volta che si si incontra un problema
alle derivate parziali.
Si consideri dunque il problema differenziale, definito come:
u (x, t) − µxx (x, t) = f (x) Ω × (0, T )
t
u(0, t) = u(L, t) = 0 t ∈ (0, T )
u(x, 0) = u (x)
0 x ∈ Ω = (0, L)
Per prima cosa è necessario passare alla formulazione debole del problema, data
da:
u ∈ V = H01 (Ω) a(u, v) = F (v) ∀v ∈ V
∑
Nh
uh (x) = φj (x)uj (t) uh (t) ∈ Vh ∀t > 0
j=1
Questa volta non si hanno dei coefficienti, ma delle vere e proprie funzioni temporali
241
uh (t) come incognite:
uj (t) : (0, T ) −→ R j = 1, 2, . . . , Nh
Inoltre, essendo la base funzione della sola coordinata spaziale l’intera derivata
temporale si scarica sull’incognita:
∂uh ∑Nh
∑
Nh ∫ L ∑
Nh ∫ L
u̇j (t) φj φi dx + uj a(φj , φi ) = f φi dx ∀i
| {z }
j=1 | 0
{z } j=1
[A]ij | {z }
0
[M]ij fi
242
non resta che procedere con la discretizzazione del dominio temporale e procedere
con un metodo risolutivo per sistemi di ODE. Ad esempio, si utilizzano i ϑ-metodi:
1 [ ]
M(uk+1 − uk ) + A ϑuk+1 + (1 − ϑ)uk = ϑf (tk+1 ) − (1 − ϑ)ftk
∆t
e ≤ chα(r) + ∆tβ(ϑ)
∆t < Cϑ h2
Vale a dire una funzione che va da un dominio noto allo spazio degli elementi finiti
che si sta considerando.
243
Tramite questa formulazione si ha che:
Dato che l’interpolante appartiene allo spazio degli elementi finiti è sempre possi-
bile considerarlo come combinazione lineare delle funzioni della base dello spazio,
cioè le lagrangiane:
∑
N +1
Πrh v(x) = v(xj )φj (x)
j=0
M r
∥u − vh ∥V ≤ c h |u|H r+1
α
244
Dato che l’interpolante Πrh u ∈ Vh ⊂ V ∩ Xhr , si può utilizzare proprio questa
funzione per ottenere un confronto:
M
∥u − uh ∥V ≤ ∥u − Πrh u∥V
α
M r
∥u − vh ∥V ≤ C h |u|H r+1
α
Teorema 12.6.2 (Errori metodi elementi finiti). Se u ∈ H p+1 (Ω) con il parametro
p > 0 si ha che:
Infine se p ≥ r:
245
Parte IV
Trasformate integrali
246
Capitolo 13
Trasformata di Fourier
+∞ [ ( ) ( )]
a0 ∑ nπ nπ
f (x) = + an cos x + bn sin x (13.1)
2 n=1
L L
247
• m=0 ∫ ∫
L L
a0
f (x)dx = dx = La0
−L −L 2
• m>0
∫ L ( ) ∫ L ( ) ( )
mπ nπ mπ
f (x) cos x dx = an cos x cos x dx = Lδn,m a0
−L L −L L L
f (x) : [0, L] −→ R
può essere sempre scritta come somma di soli seni o soli coseni. Per farlo bisogna,
prima di calcolare i coefficienti di Fourier, renderla pari, per risistemare la funzione
in modo che sia pari (per serie di coseni) oppure dispari (per la serie di seni).
1[ ]
f (x− ) + f (x+ ) x ∈ [−L, L]
2
√
cn = a2n + b2n
248
e la quantità:
nπ
ωn =
L
è detta frequenza n-esima.
In questo caso a e b non sono più coefficienti, ma vere e proprie funzioni continue,
che evolvono come:
∫ ∫
1 1
a(ω) = f (x) cos(ωx)dx b(ω) = f (x) sin(ωx)dx
π R π R
Per calcolare a(ω) non si fa altro che applicare la definizione e calcolare l’inte-
249
grale per via analitica:
∫
1
a(ω) = χ(−1,1) (x) cos(ωx)dx
π R
∫
1 1
= cos(ωx)dx
π −1
1
1
= sin(ωx)
πω −1
2
= sin(ω)
πω
Da ultimo, si ha che:
√ 2 sin(ω)
c(ω) = a2 (ω) + b2 (ω) = |a(ω)| =
π ω
250
χ(−1,1) c(ω)
2
π
1
x ω
−1 1
Inoltre si ha che:
∫ ∫
1 +∞ ( )
f (x) = f (ξ) cos ω(ξ − x) dξdω ∀x ∈ R
π 0 R
1 [ iωx ]
cos(ωx) = e + e−iωx (13.2)
2
251
Infine, sostituendo la funzione coseno con l’Eq. (13.2), si ottiene:
∫ ∫
1
f (x) = f (ξ)e−iωξ dξ eiωx dω
2π R
|R {z }
Trasformata di Fourier
252
f (x) fˆ(ω)
2
1 c
x ω
Figura 13.2: Esempio della trasformata di Fourier della funzione f (x) = e−c|x| .
Nella Fig. 13.2 è possibile avere un confronto tra la funzione originale e la sua
trasformata.
253
Teorema 13.2.2 (Proprietà di traslazione). Una traslazione viene trasformata in
una moltiplicazione per una costante opportuna. Infatti ∀ f ∈ L1 (R) e ∀ ξ0 ∈ R:
[ ] [ ]
F f (x − ξ0 ) (ω) = e−iξ0 ω F f (x) (ω)
ed è commutativo:
(f ∗ g)(x) = (g ∗ f )(x)
254
Tabella 13.1: Trasformate di Fourier delle funzioni più comuni (a, α ∈ R, β ∈ R+ )
[ ]
f (x) F f (x) (ω)
2β
e−β|x|
β 2 + ω2
1 π a|ω|
e
x2 + a2 a
sin(aω)
kχ(−a,a) (x) 2k
ω
δ0 (x) 1
1 2πδ0 (ω)
sin(αx)
πχ(−α,α) (ω)
x
√
−αx2 π ω2
e e 4α
α
255
Teorema 13.2.6 (Teorema di convoluzione). Date due funzioni f, g integrabili su
un intervallo [−L, L], in cui L < +∞, e tali per cui ∃ (f ∗ g)(x) si ha che:
[ ] [ ] [ ]
F (f ∗ g)(x) (ω) = F f (x) (ω)F g(x) (ω) = fˆ(ω)ĝ(ω)
[ ]
F −1 fˆ(ω)ĝ(ω) (x) = (f ∗ g)(x)
Teorema 13.2.7 (Trasformata della derivata). Dato n > 0 e una funzione reale
tale che f (n−1) ∈ C 0 (R), f (n−1) ∈ L1 (R) e f (n) continua a tratti su R e:
allora si ha che:
[ ]
F f (n) (x) (ω) = (iω)n fˆ(ω)
Infine si ha che:
[ ] d
F xf (x) (ω) = i fˆ(ω)
dω
Nell’ambiente MATLAB® è disponibile una funzione che permette di calcolare
la trasformata di Fourier in maniera simbolica:
F = fourier(f, t, w)
f = ifourier(F, w, t)
256
13.3 Problemi differenziali
Spesso è possibile, e sopratutto conveniente, trasformare con Fourier un problema
alle derivate parziali per trovare un metodo di risoluzione più semplice ed efficace.
Una volta trovata la soluzione è possibile tornare indietro tramite l’antitrasformata.
Purtroppo questo metodo risolutivo permette di trovare come soluzioni, le sole
funzioni trasformabili con Fourier.
|x|→+∞
u, ux −−−−−→ +∞
[ ] [ ] ∂2
F ∇2 u(x, y) (ω, y) = F uxx (x, y) + uyy (x, y) (ω, y) = −ω 2 û(ω, y) + 2 û(ω, y)
∂y
257
Mettendo tutto insieme si ottiene il problema trasformato:
û (ω, y) − ω 2 û(ω, y) = 0 (ω, y) ∈ R × (0, +∞)
yy
û(ω, 0) = fˆ(ω) ω∈R
Questo è un problema differenziale ordinario nella variabile y, molto più facile da ri-
solvere rispetto a prima. Il problema è, almeno a prima vista, non è completamente
definito: manca infatti una condizione al contorno.
In ogni caso si può trovare la soluzione generale, cioè un’espressione del tipo:
Tuttavia per rispettare l’ipotesi si deve fare in modo che la soluzione u(x, y)
sia limitata, quando y → +∞. In particolare si ha che:
û(ω, y) = k(ω)e−|ω|y
tal modo basta che il termine k(ω) sia limitato; si pongono quindi le condizioni al
contorno:
û(ω, 0) = k(ω) = fˆ(ω)
û(ω, y) = fˆ(ω)e−|ω|y
A questo punto si ha la soluzione sul dominio delle frequenze, quindi si deve so-
lo antitrasformare quanto ottenuto per concludere il problema. Si consideri la
soluzione hatu(ω, y) come prodotto di due funzioni trasformate:
û(ω, y) = fˆ(ω)ĝ(ω, y)
258
In cui la funzione:
ĝ(ω, y) = e−|ω|y
Ω = R × (0, +∞)
259
Da cui derivano i seguenti due problemi:
Dato che il dominio spaziale è infinito si hanno per forza pulsazioni ω continue
e non discrete come nei casi analizzati finora. Naturalmente anche la soluzione
u(x, t) sarà una somma continua di uω , cioè l’integrale di Fourier:
∫
u(x, t) = uω (x, t)dω
∫ R
1
a(ω) = ϕ(x) cos(ωx)dx
π R
∫
1
b(ω) = ϕ(x) sin(ωx)dx
π R
L’equazione invece: [ ] [ 2 ]
∂u ∂ u
F (ω, t) = µF (ω, t)
∂t ∂x2
A questo punto il dominio è cambiato: (ω, t) ∈ R × R+ . Il problema differenziale,
dopo la trasformazione è diventato un un problema ODE:
û(ω, t) = c(ω)e−µω
2t
260
A questo punto si applicano le condizioni iniziali:
Quindi non resta che tornare nel dominio originale: per farlo si considera la
trasformata della soluzione û(x, t) come prodotto di due trasformate:
ĝ(ω, t) = e−µω
2t
û(x, t) = ϕ̂(ω)ĝ(ω, t)
Si consideri dunque la funzione ĝ: è una trasformata notevole. Per cui si ottiene:
[ 2 ] 1 x2
g(x, t) = F −1 e−µω t (x, t) = √ e− 4µt
4πµt
1 x2
Γµ (x, t) = √ e− 4µt
4πµt
Esempio 13.3.2 (Delta di Dirac). Sia ϕ(x) = δ0 (x). Allora la soluzione è pari a:
∫
1 (x−s)2
u(x, t) = (ϕ ∗ Γµ )(x, t) = √ δ0 (s)e− 4πµt ds = Γµ (x, t)
4πµt R
261
Capitolo 14
Trasformata di Laplace
3. La trasformata F (s) è definita per valori di s tali per cui è definito l’integrale.
s = σ + iω s∈C
262
Definizione 14.0.2 (Trasformata di Laplace inversa). Data una trasformata F (s),
la sua inversa è:
∫
−1
[ ] 1 γ+iT
f (t) = L F (s) (t) = lim F (s)est ds
2πi T →+∞ γ−iT
263
f (t) F (s)
t s
f (t) F (s)
a
t s
a
Figura 14.2: Esempio della trasformata di Laplace della funzione f (x) = eat .
264
In Fig. 14.2 è possibile avere un confronto tra la funzione originale e la sua
trasformata.
∄ α, M |f (t)| ≤ M eαt
Pertanto le ipotesi del Teo. 14.1.1 non sono verificate e quindi non è possibile
trasformare la funzione.
Teorema 14.1.2 (Proprietà della linearità). Date due funzioni f, g tali che am-
mettano la trasformata di Laplace, ∀ a, b ∈ R:
[ ] [ ] [ ]
L af (t) + bg(t) (s) = aL f (t) (s) + bL g(t) (s)
[ ] [ ] [ ]
L −1 aF (s) + bG(s) (t) = aL −1 F (s) (t) + bL −1 G(s) (t)
Teorema 14.1.3 (Estremi della trasformata). Data una funzione f (t), che am-
metta la trasformata di Laplace F (s) si ha che:
Teorema 14.1.4 (Limitatezza della trasformata). Data una funzione f (t), che
ammette la trasformata di Laplace F (s) e se |f (t)| ≤ M eαt allora:
M
|F (s)| ≤ s>a
s−a
265
Tabella 14.1: Trasformate di Laplace delle funzioni più comuni (a, b, α ∈ R)
[ ]
f (t) L f (t) (s)
1
1
s
1
eat
s−a
n!
tn
sn+1
n!
tn eat
(s − a)n+1
a−b
eat − ebt
(s − a)(s − b)
α
sin(αt)
α2 + s2
s
cos(αt)
α2 + s2
α
sinh(αt)
s2 − α2
s
cosh(αt)
s2 − α2
δ0 (t − a) e−as
e−as
H(t − a)
s
∫ t
1 [ ]
f (τ )dτ L f (t) (s)
0 s
∫ +∞
1
f (t) F (σ)dσ
t s
dn [ ]
tn f (t) (−1)n F (s)
dsn
266
Teorema 14.1.5 (Primo teorema di Shift). Data una funzione f (t), che ammette
la trasformata di Laplace F (s), si ha che:
[ ]
L eat f (t) (s) = F (s − a)
[ ]
L −1 F (s − a) (t) = eat f (t)
Teorema 14.1.6 (Secondo teorema di Shift). Data una funzione f (t), che am-
mette la trasformata di Laplace F (s), si ha che:
[ ]
L H(t − a)f (t − a) (s) = e−as F (s)
[ ]
L −1 e−as F (s) (t) = H(t − a)f (t − a)
F = laplace(f, t, s)
f = ilaplace(F, s, t)
267
14.2 Trasformata di Laplace per ODE
Prima di applicare la trasformata di Laplace per risolvere dei problemi differenziali
ordinari e bene chiedersi come si comporta l’operatore quando si trasforma una
funzione derivata.
Teorema 14.2.1 (Trasformata della derivata). Sia y(t) ∈ C 0 (0, +∞) e y ′ (t)
continua a tratti su (0, k) ∀ k > 0; se:
lim y(t)e−sk = 0
k→+∞
allora:
[ ] [ ]
L y ′ (t) (s) = sL y(t) (s) − y(0) = sY (s) − y(0)
268
y ′ (t) = dy
dt
(t) ?
f (t) Sistema y(x)
L L −1
Yd = sY T (s)
F (s) Sistema Laplace Y (s) = T (s)F (s)
Figura 14.3: Schema concettuale della procedura per la risoluzione di ODE tramite l’u-
so della trasformata di Laplace. Da notare come l’operazione di derivazione si trasformi
in semplice problema algebrico.
Come si può facilmente intuire dalla formula, la risposta del sistema y(t) a un
forzamento imposto f (t) sarà data dalla convoluzione dell’antitrasformata della
funzione di trasferimento con la funzione f (t).
Cosa ben più importante, l’effetto del sistema dinamico all’interno del dominio
di Laplace, è equivalente a una semplice moltiplicazione.
La funzione di trasferimento può sempre essere riscritta come:
∏
N (s) (s − zi )
T (s) = ρ = ∏i
D(s) j (s − pi )
dove {zi } sono detti zeri della funzione e i {pi } sono chiamati poli. A seconda dei
valori degli zeri e dei poli è possibile ricavare la stabilità di un sistema dinamico;
infatti se:
ℜ{pj } < 0 ∀j
269
riassunto.
1
s2 Y (s) + 2sY (s) − 4Y (s) = + 3s + 1
s
1 3s + 1
Y (s) = + 2
s(s2 + 2s − 4) s + 2s − 4
270
Appendice A
Codici di MATLAB®
Codice A.1: Algoritmo della sostituzione in avanti per le matrici triangoli inferiori.
1 function [y] = fwsub(L, b)
2
3 % y vettore soluzione del sistema
4 % L matrice triangolare inferiore
5 % b vettore termini noti
6
7 [m , n] = size(L);
8
9 % Verifiche sui dati
10
11 flag = 0;
12 for ii = 1:n
13 for jj = (ii + 1):n
271
14 if L(ii, jj) ~= 0
15 flag = 1;
16 end
17 end
18 end
19
20 assert(prod(diag(L)) ~= 0, 'Matrice singolare')
21 assert(m == n, 'Matrice non quadrata')
22 assert(n == length(b), 'Dimensioni matrice e vettore non consistenti
')
23 assert(flag == 0, 'Matrice non triangolare inferiore')
24
25 % Algoritmo
26
27 y = zeros(n, 1);
28 y(1) = b(1)/L(1, 1);
29
30 for ii = 2:n
31 y(ii) = (b(ii) − L(ii, 1:(n − 1))*y(1:(n − 1))) / L(ii, ii);
32 end
33
34 end
Codice A.2: Algoritmo della sostituzione all’indietro per le matrici triangoli superiori.
1 function [x] = bwsub(U, y)
2
3 % x vettore soluzione del sistema
4 % U matrice triangolare superiore
5 % b vettore termini noti
6
7 [m , n] = size(U);
8
9 % Verifiche sui dati
272
10
11 flag = 0;
12 for ii = 1:n
13 for jj = 1:(ii − 1)
14 if U(ii, jj) ~= 0
15 flag = 1;
16 end
17 end
18 end
19
20 assert(prod(diag(U)) ~= 0, 'Matrice singolare')
21 assert(m == n, 'Matrice non quadrata')
22 assert(n == length(y), 'Dimensioni matrice e vettore non consistenti
')
23 assert(flag == 0, 'Matrice non triangolare superiore')
24
25 % Algoritmo
26
27 x = zeros(n, 1);
28 x(n) = y(n)/U(n, n);
29
30 for ii = (n − 1):−1:1
31 x(ii) = (y(ii) − U(ii, (ii + 1):n)*x((ii + 1):n)) / U(ii, ii);
32 end
33
34 end
273
6 % toll tollerenza sulla soluzione
7 % x soluzione del sistema
8 % niter numero delle iterazioni eseguite
9 % res residuo normalizzato
10 % incr incrementi relativi
11
12 [n, m] = size(A);
13
14 % Controllo sui dati iniziali
15
16 assert(m == n, 'Matrice non quadrata')
17 assert(all(all(A == A.')), 'Matrice non simmetrica')
18 assert(all(eig(A) > 0), 'Matrice non definta positiva')
19
20 assert(n == length(b), 'Dimensioni vettore termini noti non valido')
21 assert(n == length(x0), 'Dimensioni prima iterata non valide')
22
23 % Algoritmo
24
25 x = x0;
26 res = zeros(maxiter, 1);
27 incr = zeros(maxiter, 1);
28
29 r = b − A*x;
30
31 for ii = 1:maxiter
32 alpha = (r.'*r)/(r.'*A*r);
33 incr(ii) = norm(alpha*r)/norm(x);
34 x = x + alpha*r;
35 r = b − A*x;
36 res(ii) = norm(r)/norm(b);
37 if res(ii) < toll
38 niter = ii;
274
39 return
40 end
41 end
42
43 warning('Raggiunto numero massimo iterazioni')
44 niter = maxiter;
45
46 end
275
23 % Algoritmo
24
25 x = x0;
26 res = zeros(maxiter, 1);
27 incr = zeros(maxiter, 1);
28
29 r = b − A*x;
30 p = r;
31
32 for ii = 1:maxiter
33 alpha = (p.'*r)/(p.'*A*p);
34 incr(ii) = norm(alpha*p)/norm(x);
35 x = x + alpha*p;
36 r = b − A*x;
37 res(ii) = norm(r)/norm(b);
38 beta = (p.'*A*r)/(p.'*A*p);
39 p = r − beta*p;
40 if res(ii) < toll
41 niter = ii;
42 return
43 end
44 end
45
46 warning('Raggiunto numero massimo iterazioni')
47 niter = maxiter;
48
49 end
276
Codice A.5: Metodo di newton semplice.
1 function [xvect, it] = newton(x0, nmax, toll, fun, dfun)
2
3 % x0 punto iniziale
4 % nmax numero massimo di iterazioni
5 % toll tollerenza
6 % fun funzione da azzerare
7 % dfun derivata funzione
8 % xvect vettore contentente tutte le iterate
9 % it numero di iterazioni eseguite
10
11 xvect = zeros(nmax, 1);
12 xvect(1) = x0 − fun(x0)/dfun(x0);
13
14 for it = 2:nmax
15 xvect(it) = xvect(it − 1) − fun(xvect(it − 1))/dfun(xvect(it −
1));
16 if abs(xvect(it) − xvect(it − 1)) < toll
17 return
18 end
19 end
20
21 warning('Raggiunto numero massimo di iterazioni')
22
23 end
277
7 % dfun derivata funzione
8 % mol molteplictà zeri
9 % xvect vettore contentente tutte le iterate
10 % it numero di iterazioni eseguite
11
12 xvect = zeros(nmax, 1);
13 xvect(1) = x0 − mol*fun(x0)/dfun(x0);
14
15 for it = 2:nmax
16 xvect(it) = xvect(it − 1) − mol*fun(xvect(it − 1))/dfun(xvect(it
− 1));
17 if abs(xvect(it) − xvect(it − 1)) < toll
18 return
19 end
20 end
21
22 warning('Raggiunto numero massimo di iterazioni')
23
24 end
278
13
14 for niter = 1:nmax
15 delta = − JF(x)\F(x);
16 x = x + delta;
17 if norm(delta) < tol
18 R = delta;
19 return
20 end
21 end
22
23 warning('Raggiunto numero di iterazioni massimo')
24 R = delta;
25
26 end
279
13 x_med(ii) = mean([x(ii), x(ii + 1)]);
14 end
15
16 I = h*sum(fun(x_med));
17
18 end
280
11
12 for ii = 1:N
13 somma = somma + fun(x(ii)) + fun(x(ii + 1)) + 4*fun(mean([x(ii),
x(ii + 1)]));
14 end
15
16 I = somma*h/6;
17
18 end
281
18 end
282
9 t = (a:h:b).'; % Nodi temporali
10 u = zeros(size(t)); % Azzeramento vettore soluzione
11
12 u(1) = y_0;
13
14 for n = 1:(length(t) − 1)
15 fun = @(w) − w + u(n) + 0.5*h.*(f(w, t(n + 1)) + f(u(n), t(n)));
16 u(n + 1) = fsolve(fun, u(n));
17 clc
18 end
19
20 end
283