Dispense del corso di Laboratorio di Metodi Numerici per le

Dispense del corso
Laboratorio di Metodi Numerici per le
Equazioni Differenziali
Dott. Marco Caliari
a.a. 2009/10
aggiornate al 6 ottobre 2010
Questi appunti non hanno nessuna pretesa di completezza. Sono solo alcune note ed esercizi che affiancano il corso di Metodi Numerici per le Equazioni Differenziali. Sono inoltre da considerarsi in perenne “under revision”
e pertanto possono contenere discrepanze, inesattezze o errori.
Indice
0 Preliminari
0.1 Metodi iterativi per sistemi di equazioni lineari . . . .
0.2 Metodi di Richardson . . . . . . . . . . . . . . . . . . .
0.2.1 Metodo del gradiente precondizionato . . . . . .
0.2.2 Metodo del gradiente coniugato precondizionato
0.2.3 Test d’arresto . . . . . . . . . . . . . . . . . . .
0.3 Memorizzazione di matrici sparse . . . . . . . . . . . .
0.3.1 Alcuni comandi per matrici sparse . . . . . . . .
0.4 Metodo di Newton per sistemi di equazioni non lineari
0.4.1 Metodo di Newton modificato . . . . . . . . . .
0.5 Esponenziale di matrice . . . . . . . . . . . . . . . . .
0.5.1 Formula delle variazioni delle costanti . . . . .
0.5.2 Calcolo di exp(A) . . . . . . . . . . . . . . . . .
0.6 Esercizi . . . . . . . . . . . . . . . . . . . . . . . . . .
I
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
ODEs(Equazioni differenziali ordinarie)
6
6
6
7
8
9
10
10
11
12
12
12
13
16
18
1 theta-metodo
1.1 Caso lineare . . . . . . . . . . . . . . . . . . . . . . . . . . . .
1.2 Equazioni di ordine superiore al primo . . . . . . . . . . . . .
1.3 Verifica della correttezza dell’implementazione . . . . . . . . .
19
20
20
21
2 Metodi di Runge–Kutta
2.1 Implementazione dei metodi di Runge–Kutta
2.2 Metodi di Runge–Kutta embedded . . . . . .
2.2.1 Passo di integrazione variabile . . . .
2.2.2 Strutture in GNU Octave . . . . . .
24
24
24
26
26
3 Integratori esponenziali
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
27
3
4
INDICE
4 Esercizi
29
II
32
BVPs(Problemi ai limiti)
5 Metodo di shooting
34
5.1 Metodo di bisezione . . . . . . . . . . . . . . . . . . . . . . . . 34
5.2 Metodo di Newton . . . . . . . . . . . . . . . . . . . . . . . . 35
5.3 Problema ai limiti con frontiera libera . . . . . . . . . . . . . . 36
6 Differenze finite
6.1 Costruzione delle matrici delle differenze finite
6.2 Condizioni di Dirichlet . . . . . . . . . . . . .
6.3 Condizioni di Neumann . . . . . . . . . . . . .
6.4 Differenze finite centrate del secondo ordine .
6.5 Convergenza per un problema modello . . . .
6.5.1 Unicità . . . . . . . . . . . . . . . . . .
6.5.2 Esistenza . . . . . . . . . . . . . . . .
6.5.3 Regolarità . . . . . . . . . . . . . . . .
6.5.4 Consistenza . . . . . . . . . . . . . . .
6.5.5 Esistenza ed unicità . . . . . . . . . . .
6.5.6 Proprietà di Ah . . . . . . . . . . . . .
6.5.7 Stabilità . . . . . . . . . . . . . . . . .
6.5.8 Convergenza . . . . . . . . . . . . . . .
6.5.9 Un esempio: l’equazione della catenaria
6.6 Norme . . . . . . . . . . . . . . . . . . . . . .
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
38
38
39
40
40
41
41
42
42
43
44
44
44
45
45
46
7 Equazione di Poisson
47
7.1 Equazione di Poisson bidimensionale . . . . . . . . . . . . . . 47
7.1.1 Condizioni al bordo di Dirichlet . . . . . . . . . . . . . 47
7.1.2 Condizioni al bordo miste . . . . . . . . . . . . . . . . 49
8 Metodi variazionali
8.1 Formulazione variazionale di un problema modello
8.1.1 Metodo di approssimazione variazionale . .
8.1.2 Estensione al caso bidimensionale . . . . .
8.2 Metodi spettrali . . . . . . . . . . . . . . . . . . .
8.2.1 Trasformata di Fourier . . . . . . . . . . .
8.2.2 Trasformata di Fourier discreta . . . . . .
8.3 Metodi di collocazione . . . . . . . . . . . . . . .
8.3.1 Condizioni al bordo . . . . . . . . . . . . .
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
51
51
53
57
57
58
59
64
65
5
INDICE
9 Esercizi
68
III
PDEs(Equazioni alle derivate parziali)
70
10 Equazione del calore
10.1 Equazione del calore con dati iniziali e condizioni ai limiti
10.1.1 Esistenza di una soluzione . . . . . . . . . . . . . .
10.1.2 Unicità della soluzione . . . . . . . . . . . . . . . .
10.2 Metodo delle linee . . . . . . . . . . . . . . . . . . . . . . .
10.2.1 Differenze finite . . . . . . . . . . . . . . . . . . . .
10.2.2 Elementi finiti . . . . . . . . . . . . . . . . . . . . .
10.3 Esercizi . . . . . . . . . . . . . . . . . . . . . . . . . . . .
.
.
.
.
.
.
.
.
.
.
.
.
.
.
71
71
71
73
74
74
75
77
11 Equazione del trasporto
11.1 Linee caratteristiche . . . . . . . . . . .
11.2 Differenze finite . . . . . . . . . . . . . .
11.2.1 Eulero esplicito/upwind . . . . .
11.2.2 Condizione CFL . . . . . . . . . .
11.2.3 Equazioni modificate . . . . . . .
11.2.4 Formulazione di flusso . . . . . .
11.2.5 Verifica dell’ordine . . . . . . . .
11.2.6 Estensione . . . . . . . . . . . . .
11.2.7 Condizioni al bordo . . . . . . . .
11.2.8 Problemi di convezione-diffusione
11.3 Esercizi . . . . . . . . . . . . . . . . . .
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
78
78
80
80
81
83
84
85
85
85
85
86
Bibliografia
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
88
Capitolo 0
Preliminari
0.1
Metodi iterativi per sistemi di equazioni
lineari
I metodi iterativi per la soluzione del sistema lineare
Ax = b
(1)
si basano sull’idea di calcolare la soluzione come limite di una successione di
vettori
x = lim x(l) .
l→∞
Una strategia generale per costruire la successione {x(l) }l è basata sullo splitting A = P − M , ove P è non singolare. Assegnato x(1) , il termine x(l+1) è
calcolato ricorsivamente come
P x(l+1) = M x(l) + b,
l≥1
(2)
Posto e(l) = x − x(l) , si ha
e(l) = Bel−1 ,
B = P −1 M = I − P −1 A ,
ove B è chiamata matrice di iterazione.
Lemma 1. Si ha liml→∞ e(l) = 0 per ogni e(1) se e solo se liml→∞ B l = 0,
cioè se e solo se ρ(B) < 1.
0.2
Metodi di Richardson
Indicato con r(l) il residuo
r(l) = b − Ax(l) = Ax − Ax(l) = A(x − x(l) ) = Ae(l) ,
6
7
0.2. METODI DI RICHARDSON
il metodo (2) può essere riscritto come
P (x(l+1) − x(l) ) = r(l) .
(3)
In questo contesto, P viene chiamata matrice di precondizionamento o precondizionatore di A e viene scelta in modo che la matrice di iterazione
B = I − P −1 A abbia un raggio spettrale minore di 1 e la risoluzione di
(3) sia “facile”.
Una generalizzazione dello schema (3) è il metodo di Richardson: dato
x(1) , x(l+1) è calcolato ricorsivamente come
P (x(l+1) − x(l) ) = αr(l) ,
ove α è un opportuno parametro di accelerazione. Dati x(1) e r(1) = b−Ax(1) ,
l’algoritmo per calcolare x(l+1) è
P z (l) = r(l)
x(l+1) = x(l) + αz (l)
(4)
r(l+1) = r(l) − αAz (l)
Il costo di un’iterazione è dato essenzialmente dalla risoluzione di un sistema
lineare P z (l) = r(l) facile e dal prodotto matrice-vettore Az (l) . Tali metodi
risulteranno particolarmente vantaggiosi per matrici sparse, in cui il numero
di elementi diversi da zero è O(N ) piuttosto che O(N 2 ) (e dunque il costo
di un prodotto matrice-vettore è O(N )), se l’ordine della matrice è N .
Il calcolo del residuo r(l+1) = r(l) − αAz (l) (invece di r(l+1) = b − Ax(l+1) )
permette di ridurre la propagazione, attraverso il prodotto matrice-vettore,
degli errori, in quanto il vettore z (l) , contrariamente a x(l+1) , diminuisce in
modulo al crescere di l.
0.2.1
Metodo del gradiente precondizionato
Siano A e P simmetriche e definite positive. Il metodo di Richardson può
essere generalizzato con una scelta dinamica del parametro di accelerazione,
prendendo α = αl in modo tale che
p
kx − x(l+1) kA , kykA = y T Ay
sia minima. Si ha
kx − x(l+1) k2A = (x − x(l) − αl z (l) )T A(x − x(l) − αl z (l) ) =
T
T
= αl2 z (l) Az (l) − 2αl z (l) A(x − x(l) ) + (x − x(l) )T A(x − x(l) )
8
CAPITOLO 0. PRELIMINARI
e dunque il minimo è dato dalla scelta
T
αl =
z (l) r(l)
T
z (l) Az (l)
.
Il metodo ottenuto si chiama metodo del gradiente precondizionato. Dati x(1)
e r(1) , l’algoritmo per calcolare x(l+1) è
P z (l) = r(l)
T
z (l) r(l)
αl =
T
z (l) Az (l)
(l+1)
x
= x(l) + αl z (l)
(5)
r(l+1) = r(l) − αl Az (l)
Nel caso si scelga P = I, si ottiene il metodo del gradiente (noto anche come
steepest descent).
0.2.2
Metodo del gradiente coniugato precondizionato
Siano A e P simmetriche e definite positive. Il metodo del gradiente coniugato
precondizionato è una generalizzazione del metodo di Richardson in cui
x(l+1) = x(l) + αl p(l)
ove i {p(l) }l sono coniugati, cioè soddisfano
T
p(i) Ap(j) = 0,
i 6= j
Per soddisfare questa proprietà è necessaria l’introduzione di un ulteriore
parametro βl . Dati x(1) , r(1) , P z (1) = r(1) e p(1) = z (1) , l’algoritmo per
calcolare x(l+1) è
T
z (l) r(l)
αl =
T
p(l) Ap(l)
x(l+1) = x(l) + αl p(l)
r(l+1) = r(l) − αl Ap(l)
P z (l+1) = r(l+1)
T
βl+1 =
p(l+1)
z (l+1) r(l+1)
T
z (l) r(l)
(l+1)
=z
+ βl+1 p(l)
(6)
0.2. METODI DI RICHARDSON
9
Teorema 1. Il metodo del gradiente coniugato applicato ad una matrice di
ordine N converge in al più N iterazioni (in aritmetica esatta).
Dimostrazione. La dimostrazione (omessa) si basa essenzialmente sul fatto
che p(1) , . . . , p(N ) sono vettori linearmente indipendenti e non ce ne possono
essere più di N .
Per questo motivo, tale metodo è detto semiiterativo.
Stima dell’errore
Vale la seguente stima dell’errore:
!l−1
Ãp
−1 A) − 1
cond
(P
2
ke(l) kA ≤ 2 p
ke(1) kA
−1
cond2 (P A) + 1
dalle quale si osserva che
• la stima d’errore decresce in ogni caso, poiché il numeratore è più
piccolo del denominatore;
• in particolare, nel caso P = I;
• tanto più è piccolo il numero di condizionamento di P −1 A, tanto più il
metodo ha convergenza veloce;
• nel caso limite di P = A, si ha ke(l) kA ≤ 0.
0.2.3
Test d’arresto
Un primo stimatore è costituito dal residuo: si arresta cioè il metodo iterativo
quando
kr(l) k ≤ tol · kbk .
Infatti, dalla precedente si ricava
ke(l) k
≤ tol · cond(A) .
kxk
Una modifica consiste in
kr(l) k ≤ tol · kr(1) k ,
(7)
che coincide con il precedente nel caso in cui come x(1) venga scelto il vettore
di zeri.
10
CAPITOLO 0. PRELIMINARI
0.3
Memorizzazione di matrici sparse
Sia A una matrice sparsa di ordine N con m elementi diversi da zero. Esistono
molti formati di memorizzazione di matrici sparse. Quello usato da GNU
Octave è il Compressed Column Storage (CCS). Consiste di tre array: un
primo, data, di lunghezza m contenente gli elementi diversi da zero della
matrice, ordinati prima per colonna e poi per riga; un secondo, ridx, di
lunghezza m contenente gli indici di riga degli elementi di data; ed un terzo,
cidx, di lunghezza N +1, il cui primo elemento è 0 e l’elemento i+1-esimo è il
numero totale di elementi diversi da zero nelle prime i colonne della matrice.
Per esempio, alla matrice


1 0 0 0
0 2 3 0

A=
4 0 5 6
0 0 0 7
corrispondono i vettori
data = [1, 4, 2, 3, 5, 6, 7]
ridx = [1, 3, 2, 2, 3, 3, 4]
cidx = [0, 2, 3, 5, 7]
In GNU Octave, il formato CCS e l’implementazione del prodotto matrice-vettore sono automaticamente usati dalla function sparse e dall’operatore *, rispettivamente.
0.3.1
Alcuni comandi per matrici sparse
• Il comando speye(N) genera la matrice identità di ordine N .
• Il comando spdiags(v,0,N,N), ove v è un vettore colonna, genera la
matrice diagonale di ordine n avente v in diagonale. Se la dimensione
di v è minore di n, la diagonale viene riempita con zeri posti dopo il
vettore v. Se invece la dimensione di v è maggiore di N , vengono usate
solo le prime N componenti di v.
Sia V la matrice

v11
 v21

V =  ..
 .
v12
v22
..
.

v13
v23 

.. 
. 
vN 1 vN 2 vN 3
0.4. METODO DI NEWTON PER SISTEMI DI EQUAZIONI NON LINEARI11
Il comando spdiags(V,-1:1,N,N) genera la matrice
v12
v11


0
 .
 ..

0
0

0.4
v23 0
v22 v33
0
0
...
...
...
...
...
...
0
0
..
.
..
.








vN 3 
v21 v32
... ...
. . . 0 vN −2 1 vN −1 2
... ...
0
vN −1 1 vN 2
Metodo di Newton per sistemi di equazioni non lineari
Consideriamo il sistema di equazioni non lineari

f1 (x1 , x2 , . . . , xN ) = 0




 f2 (x1 , x2 , . . . , xN ) = 0
..


.



fN (x1 , x2 , . . . , xN ) = 0
che può essere riscritto, in forma compatta,
f (x) = 0 .
Dato x(1) , il metodo di Newton per calcolare x(l+1) è
J (l) δx(l) = −f (x(l) )
x(l+1) = x(l) + δx(l)
(8)
ove J (l) è la matrice Jacobiana, definita da
(l)
Jij =
∂fi (x(l) )
(l)
∂xj
Il criterio d’arresto solitamente usato è
kδx(l) k ≤ tol .
.
(9)
12
CAPITOLO 0. PRELIMINARI
0.4.1
Metodo di Newton modificato (o inesatto)
Il metodo di Newton (8) richiede il calcolo della matrice Jacobiana e la sua
“inversione” ad ogni passo k. Questo potrebbe essere troppo oneroso. Una
strategia per ridurre il costo computazionale è usare sempre la stessa matrice Jacobiana J (1) , oppure aggiornarla solo dopo un certo numero di iterazioni. In tal modo, per esempio, è possibile usare la stessa fattorizzazione
L(l) U (l) per più iterazioni successive. Un’altra strategia potrebbe essere quella di usare una matrice Jacobiana più semplice (per esempio quella risultante
dall’eventuale parte lineare della funzione).
0.5
Esponenziale di matrice
Data una matrice quadrata A ∈ RN ×N , si definisce
exp(A) =
∞
X
Aj
j=0
j!
.
Tale serie converge per qualunque matrice A, essendo A un operatore lineare
tra spazi di Banach e avendo la serie esponenziale raggio di convergenza ∞.
Se A e B sono permutabili (cioè AB = BA), allora
exp(A + B) = exp(A) exp(B) .
0.5.1
Formula delle variazioni delle costanti
Data l’equazione differenziale
(
y ′ (t) = ay(t) + b(y(t)),
y(t0 ) = y0
t>0
(10)
y(t) ∈ R, la soluzione può essere scritta analiticamente mediante la formula
delle variazioni delle costanti
Z t
(t−t0 )a
e(t−τ )a b(y(τ ))dτ
(11)
y(t) = e
y0 +
t0
Infatti, si ha
′
(t−t0 )a
y (t) = ae
y0 + a
Z
t
t0
e(t−τ )a b(y(τ ))dτ + e(t−t)a b(y(t)) = ay(t) + b(y(t))
13
0.5. ESPONENZIALE DI MATRICE
Si osservi che
Z
t
(t−τ )a
e
t0
¯t
1 (t−τ )a ¯¯
−ae
dτ = − e
¯ =
a
t0
t0
(t−t0 )a
¡
¢
1
e
−1
= − 1 − e(t−t0 )a = (t − t0 )
=
a
(t − t0 )a
= (t − t0 )ϕ1 ((t − t0 )a) ,
1
dτ = −
a
Z
t
(t−τ )a
ove
∞
ez − 1 X z j
=
.
ϕ1 (z) =
z
(j
+
1)!
j=0
(12)
e, analogamente,
Z
t
t0
e(t−τ )a (τ − t0 )dτ = (t − t0 )2 ϕ2 ((t − t0 )a)
ove
∞
ez − 1 − z X z j
ϕ2 (z) =
=
.
z2
(j + 2)!
j=0
Consideriamo ora un sistema differenziale
(
y ′ (t) = Ay(t) + b(y(t)),
y(t0 ) = y 0
t>0
(13)
(14)
Ancora, la soluzione esplicita può essere scritta come
y(t) = exp((t − t0 )A)y 0 +
0.5.2
Z
t
t0
exp((t − τ )A)b(y(τ ))dτ
Calcolo di exp(A)
Come per la risoluzione di sistemi lineari, non esiste il modo per calcolare
exp(A), ma diversi modi, ognuno adatto a particolari situazioni.
Matrici piene, di modeste dimensioni
Questi metodi si applicano, in pratica, a quelle matrici per le quali si usano
i metodi diretti per la risoluzione di sistemi lineari.
14
CAPITOLO 0. PRELIMINARI
Decomposizione spettrale Se la matrice è diagonalizzabile, cioè A =
V DV −1 , allora exp(A) = V exp(D)V −1 , ove exp(D) è la matrice diagonale
con elementi ed1 , ed2 , . . . , edN . Basta infatti osservare che
A2 = (V DV −1 )2 = (V DV −1 )(V DV −1 ) = V D2 V −1
e scrivere exp(A) come serie di Taylor. La decomposizione spettrale di una
matrice costa, in generale, O(N 3 ). Si ottiene in GNU Octave con il comando
eig.
Approssimazione razionale di Padé Si considera un’approssimazione
razionale della funzione esponenziale
a1 z p−1 + a2 z p−2 + . . . + ap
e ≈
,
b1 z q−1 + b2 z q−2 + . . . + bq
z
(15)
ove bq = 1 per convenzione. Essa è chiamata diagonale quando p = q. Si
può dimostrare che le approssimazioni diagonali sono le più efficienti. Fissato
il grado di approssimazione, si sviluppa in serie di Taylor la funzione esponenziale e si fanno coincidere quanti più coefficienti possibile. Per esempio,
fissiamo p = q = 2. Si ha allora
¶
µ
z2 z3
+
+ . . . (b1 z + 1) = a1 z + a2
1+z+
2
6
z2
b1 z + 1 + b1 z 2 + z +
+ o(z 2 ) = a1 z + a2
2
da cui




1 = a2
b 1 + 1 = a1


 b1 + 1 = 0
2
L’approssimazione di Padé si estende banalmente al caso matriciale. Considerando sempre il caso p = q = 2, si ha
exp(A) ≈ B = (b1 A + I)−1 (a1 A + a2 I) ,
cioè B è soluzione del sistema lineare (b1 A + I)B = a1 A + a2 I.
L’approssimazione di Padé è accurata solo quando |z| < 1/2 (o, nel caso
matriciale, kAk2 < 1/2). Per la funzione esponenziale esiste una tecnica,
chiamata scaling and squaring che permette di aggirare il problema. Si usa
infatti la proprietà
´j
¡
¢2 ³
j 2
ez = ez/2 = ez/2
.
15
0.5. ESPONENZIALE DI MATRICE
Se |z| > 1/2, allora |z|/2j < 1/2 per j > log2 (|z|) + 1. Si calcola dunque
j
l’approssimazione di Padé di ez/2 e poi si eleva al quadrato j volte. Per la
funzione ϕ1 vale
³z ´
1
ϕ1 (z) = (ez/2 + 1)ϕ1
2
2
Anche l’approssimazione di Padé matriciale ha costo O(N 3 ). In GNU Octave
si usa una variante di questa tecnica nel comando expm.
Matrici sparse, di grandi dimensioni
I metodi visti nel paragrafo precedente ignorano l’eventuale sparsità delle
matrici. Inoltre, negli integratori esponenziali, non è mai richiesto di calcolare
esplicitamente funzioni di matrice, ma solo funzioni di matrice applicate a
vettori, cioè exp(A)v (è l’analoga differenza tra calcolare A−1 e A−1 v). Si
possono allora usare dei metodi iterativi.
Metodo di Krylov Mediante la tecnica di Arnoldi è possibile, tramite
prodotti matrice-vettore, decomporre A in A ≈ Vm Hm VmT , ove Vm ∈ Rn×m ,
VmT Vm = I, Vm e1 = v e Hm è matrice di Hessenberg di ordine m (con m ≪
N ). Allora AVm ≈ Vm Hm e
exp(A)v ≈ Vm exp(Hm )e1 .
Il calcolo di exp(Hm ) è fatto mediante l’approssimazione di Padé. Il costo
della tecnica di Arnoldi è O(N m2 ) se A è matrice sparsa. È necessario inoltre
memorizzare la matrice Vm .
Interpolazione su nodi di Leja Se il polinomio pm (z) interpola ez nei
nodi ξ0 , ξ1 , . . . , ξm , allora pm (A)v è una approssimazione di exp(A)v. È una
buona approssimazione se i nodi sono buoni (non equispaziati, per esempio)
e se sono contenuti nell’involucro convesso dello spettro di A. È difficile
stimare a priori il grado di interpolazione m necessario. È conveniente usare
la formula di interpolazione di Newton
pm−1 (z) = d1 + d2 (z − ξ1 ) + d3 (z − ξ1 )(z − ξ2 ) + · . . . · +dm (z − ξ1 ) · · · (z − ξm−1 )
ove {di }i sono le differenze divise. Tale formula si può scrivere, nel caso
matriciale,
!
Ãm−1
Y
(A − ξi I) v = (A − ξm−1 )wm−1
pm−1 (A)v = pm−2 v + dm wm , wm =
i=1
16
CAPITOLO 0. PRELIMINARI
Dunque, la complessità è O(N m) è richiesta la memorizzazione di un solo
vettore w.
Quali nodi usare? I nodi di Chebyshev, molto buoni per l’interpolazione,
non possono essere usati, in quanto non permettono un uso efficiente della
formula di interpolazione di Newton (cambiano tutti al cambiare del grado).
I nodi di Leja sono distribuiti asintoticamente come i nodi di Chebyshev e,
dati i primi m − 1, ξm è il nodo per cui
m−1
Y
i=1
|ξm − ξi | = max
ξ∈[a,b]
m−1
Y
i=1
|ξ − ξi | ,
ove l’intervallo [a, b] è in relazione con lo spettro di A, per esempio [a, b] =
σ(A) ∩ {y = 0}. Il primo nodo coincide, solitamente, con l’estremo dell’intervallo [a, b] di modulo massimo. È chiaro che l’insieme dei primi m nodi
di Leja coincide con l’unione di {ξm } con l’insieme dei primi m − 1 nodi di
Leja.
0.6
Esercizi
1. Implemetare le functions [data,ridx,cidx] = full2ccs(A) e [A] =
ccs2full(data,ridx,cidx) e le functions che, dati data, ridx e cidx,
implementano i prodotti matrice vettore Ax e AT x.
2. Si risolvano 6 sistemi lineari con le matrici di Hilbert di ordine N =
4, 6, 8, 10, 12, 14 (hilb(N)) e termine noto scelto in modo che la soluzione esatta sia il vettore [1, 1, . . . , 1]T usando il comando \ di GNU
Octave, il metodo del gradiente precondizionato e il metodo del gradiente coniugato precondizionato. Per questi ultimi due, si usi una
tolleranza pari a 10−6 , un numero massimo di iterazioni pari a 2000,
il precondizionatore diagonale e un vettore iniziale x(1) di zeri. Si riporti, per ogni N , il numero di condizionamento della matrice, l’errore
in norma infinito rispetto alla soluzione esatta e il numero di iterazioni
dei metodi iterativi.
3. Risolvere il sistema non lineare
(
f1 (x1 , x2 ) = x21 + x22 − 1 = 0
f2 (x1 , x2 ) = sin(πx1 /2) + x32 = 0
con il metodo di Newton (8). Si usi una tolleranza pari a 10−6 , un
numero massimo di iterazioni pari a 150 e un vettore iniziale x(1) =
0.6. ESERCIZI
17
[1, 1]T . Si risolva lo stesso sistema non lineare usando sempre la matrice
Jacobiana relativa al primo passo e aggiornando la matrice Jacobiana
ogni r iterazioni, ove r è il più piccolo numero di iterazioni che permette
di ottenere la stessa soluzione con la tolleranza richiesta calcolando solo
due volte la matrice Jacobiana.
4. Si implementi una function [a,b] = padeexp(p) che restituisce i coefficienti dell’approssimazione razionale di Padé (15) (con p = q) per la
funzione esponenziale.
Parte I
ODEs
(Equazioni differenziali
ordinarie)
18
Capitolo 1
θ-metodo
Consideriamo il sistema di ODEs
 ′
y1 (t) = f1 (t, y1 (t), y2 (t), . . . , yd (t))




 y2′ (t) = f2 (t, y1 (t), y2 (t), . . . , yd (t))
..


.


 ′
yd (t) = fd (t, y1 (t), y2 (t), . . . , yd (t))
con dato iniziale

y1 (t0 ) = y1 0




 y2 (t0 ) = y2 0
..


.



yd (t0 ) = yd 0
che può essere riscritto, in forma compatta,
(
y ′ (t) = f (t, y(t))
y(t0 ) = y 0
(1.1)
Notiamo come il sistema non autonomo (1.1) può essere ricondotto ad uno
autonomo
 ′
y0 (t) = 1



 y ′ (t) = f (y (t), y(t))
0
(1.2)

y
(t
)
=
t
0 0
0



y(t0 ) = y 0
ponendo y0 (t) = t.
Il θ-metodo per il sistema (1.2) si scrive
y n+1 − y n
= θf (y n ) + (1 − θ)f (y n+1 )
h
19
(1.3)
20
CAPITOLO 1. THETA-METODO
ove tn+1 = tn + h e y n ≈ y(tn ). Chiaramente, il θ metodo si riduce al metodo
di Eulero esplicito per θ = 1, al metodo di Eulero implicito per θ = 0 e
al metodo di Crank–Nicolson per θ = 1/2. Nel caso implicito (θ 6= 1), ad
ogni passo si deve risolvere un sistema di equazioni non lineari F (x) = 0,
x = y n+1 , ove
F (x) = x − h(1 − θ)f (x) − y n − hθf (y n ).
La matrice Jacobiana associata (utile per l’applicazione del metodo di Newton) è
∂f (x)
Jij (x) = I − h(1 − θ) i
.
∂xj
Ovviamente, si può scegliere un metodi di Newton modificato. Come vettore
iniziale per il calcolo di y n+1 si può prendere la soluzione al passo precedente
yn.
1.1
Caso lineare
Il caso più frequente è quello lineare autonomo
(
y ′ (t) = Ay(t)
y(t0 ) = y 0
con passo di integrazione h costante. In tal caso, il metodo si scrive
(I − h(1 − θ)A)y n+1 = y n + hθAy n
Nel caso implicito, si tratta dunque di risolvere un sistema lineare di matrice
I − h(1 − θ)A ad ogni passo. Pertanto, per problemi di piccola dimensione,
è conveniente precalcolare la fattorizzazione LU della matrice. Altrimenti, si
può considerare un metodo iterativo, ove si scelga come vettore iniziale per
il calcolo di y n+1 la soluzione al passo precedente y n .
1.2
Equazioni di ordine superiore al primo
Le equazioni differenziali di ordine d del tipo


y (d) (t) = f (t, y(t), y ′ (t), . . . , y (d−1) (t))




y(t ) = y0,0


 ′ 0
y (t0 ) = y0,1


..


.



 (d−1)
y
(t0 ) = y0,d−1
1.3. VERIFICA DELLA CORRETTEZZA DELL’IMPLEMENTAZIONE21
si possono ricondurre ad un sistema di ODEs di ordine uno, mediante la
sostituzione

y1 (t) = y(t)




 y2 (t) = y ′ (t)
..

.




yd−1 (t) = y (d−1) (t)
dando cosı̀ luogo a
(
y ′ (t) = f (t, y(t))
y(t0 ) = [y0,0 , y0,1 , . . . , y0,d−1 ]T
ove
f (t, y(t)) = [y2 (t), y3 (t), . . . , yd (t), f (t, y1 (t), y2 (t), . . . , yd−1 (t))]T
1.3
Verifica della correttezza dell’implementazione
Supponiamo di aver implementato un metodo di ordine p per la soluzione del
sistema differenziale
(
y ′ (t) = f (y(t))
y(t0 ) = y 0
e di volerne testare la corretta implementazione. L’idea è quella di creare
una soluzione artificiale x(t), inserirla nell’equazione e calcolarne il residuo
x′ (t) − f (x(t)) = g(t)
A questo punto, si risolve il sistema differenziale
(
y ′ (t) = f (y(t)) + g(t) = f̂ (t, y(t))
y(t0 ) = x(t0 )
fino ad un tempo t0 + t∗ fissato, con due discretizzazioni di passo costante
h1 = t∗ /m1 e h2 = t∗ /m2 , rispettivamente. Si avranno errori finali em1 ,h1 =
ky m1 ,h1 − x(t0 + t∗ )k = Chp1 e em2 ,h2 = ky m2 ,h2 − x(t0 + t∗ )k = Chp2 . Si ha
dunque
µ ¶p
em2 ,h2
h2
=
,
em1 ,h1
h1
22
CAPITOLO 1. THETA-METODO
da cui
log em2 ,h2 − log em1 ,h1 = p(log h2 − log h1 ) = −p(log m2 − log m1 ) .
Dunque, rappresentando in un grafico logaritmico-logaritmico l’errore in dipendenza dal numero di passi, la pendenza della retta corrisponde all’ordine
del metodo, cambiato di segno. Tale verifica è valida anche nel caso di passi
non costanti.
Nel caso f (y(t)) sia particolarmente complicato, invece di calcolare il
residuo, si può calcolare una soluzione di riferimento y m̄,h̄ e poi confrontare
con essa le soluzioni y m1 ,h1 e y m2 ,h2 , ove m1 , m2 ≪ m̄. In questo caso,
però, si può mostrare solo che il metodo converge con l’ordine giusto ad una
soluzione, non necessariamente quella giusta.
Falsa superconvergenza
Supponiamo che la soluzione esatta di un problema differenziale sia y = 1 e
che, discretizzandolo con m = 1, 2, 4, 8 passi si ottengano le soluzioni y1 = 17,
y2 = 9, y4 = 5 e y8 = 3. Supponiamo poi di averne calcolato una soluzione
di riferimento con m = 16 passi ȳ = y16 = 2. In Figura 1.1 il grafico degli
errori che si ottiene.
1.3. VERIFICA DELLA CORRETTEZZA DELL’IMPLEMENTAZIONE23
100
rispetto alla soluzione esatta
rispetto alla soluzione di riferimento
errore
10
1
0.1
1
2
4
m
Figura 1.1: Ordine dei metodi
8
Capitolo 2
Metodi di Runge–Kutta
2.1
Implementazione dei metodi di Runge–
Kutta
Consideriamo un metodo Runge–Kutta esplicito

i−1
X




ξ
=
y
+
h
ai,j f (tn + cj h, ξ j ),
n

 i
i = 1, . . . , ν
j=1
ν

X



y
=
y
+
h
bj f (tn + cj h, ξ j )

n
 n+1
j=1
Ai fini dell’implementazione, per evitare di calcolare più volte la funzione f
negli stessi punti, si usa lo schema

i−1
X




k = f (tn + ci h, y n + h
ai,j kj ), i = 1, . . . , ν

 i
j=1
2.2




y
= yn + h

 n+1
ν
X
bj k j
j=1
Metodi di Runge–Kutta embedded
Supponiamo di avere un metodo di Runge–Kutta esplicito di ordine p − 1
il cui tableaux è riportato nella Tabella 2.1 e un altro metodo di Runge–
Kutta di ordine p il cui tableaux è riportato nella Tabella 2.2. È chiaro che,
dopo aver costruito il primo metodo, con una sola nuova valutazione della
funzione f si può costruire il secondo metodo. Una tale coppia di metodi si
24
2.2. METODI DI RUNGE–KUTTA EMBEDDED
0
c2
c3
..
.
cν−1
a2,1
a3,1
..
.
a3,2
..
.
...
aν−1,1
b1
aν−1,2
b2
...
...
25
aν−1,ν−2
bν−2
bν−1
Tabella 2.1: Metodo di Runge–Kutta di ordine p − 1.
0
c2
c3
..
.
a2,1
a3,1
..
.
a3,2
..
.
...
cν−1
cν
aν−1,1
aν,1
aν−1,2
aν,2
...
...
aν−1,ν−2
aν,ν−2
aν,ν−1
b̂1
b̂2
...
b̂ν−2
b̂ν−1
b̂ν
Tabella 2.2: Metodo di Runge–Kutta di ordine p.
dice embedded e si scrive di solito un unico tableaux, come nella Tabella 2.3.
Consideriamo il sistema differenziale
0
c2
c3
..
.
cν
a2,1
a3,1
..
.
a32
..
.
...
aν,1
b1
aν,2
b2
...
...
aν,ν−1
bν−1
b̂1
b̂2
...
b̂ν−1
b̂ν
Tabella 2.3: Metodi di Runge–Kutta embedded di ordine p − 1 e p.
(
ỹ ′ (t) = f (t, ỹ(t))
ỹ(tn ) = y (p)
n
(p)
ove y n è l’approssimazione di y(tn ) ottenuta con il metodo Runge–Kutta di
ordine p. Si ha allora
(p)
(p−1)
(p)
(p−1)
ky n+1 − y n+1 k = ky n+1 − ỹ(tn+1 ) + ỹ(tn+1 ) − y n+1 k / Chpn+1 ,
(2.1)
26
CAPITOLO 2. METODI DI RUNGE–KUTTA
ove hn+1 = tn+1 − tn è il passo di integrazione e Chp−1
n+1 è l’errore di troncamento locale del metodo di ordine p − 1. Se si vuole controllare tale errore
si può allora imporre, ad ogni passo, che
(p)
(p−1)
(p)
ky n+1 − y n+1 k ≤ tola + ky n+1 ktolr ,
(p)
(p−1)
rifiutando y n+1 (e y n+1 ) nel caso la disuguaglianza non sia soddisfatta e
calcolando un nuovo passo di integrazione minore del precedente.
2.2.1
Passo di integrazione variabile
Per calcolare il successivo passo di integrazione hn+2 , imponiamo che valga
(p)
Chpn+2 = tola + ky n+1 ktolr
e, ricavando 1/C da (2.1), si ha
hn+2 =
Ã
(p)
tola + ky n+1 ktolr
(p)
(p−1)
ky n+1 − y n+1 k
!1/p
· hn+1 .
Per evitare che il passo di integrazione cambi troppo bruscamente, si può
adottare una correzione del tipo


Ã
!1/p 
(p)
tol
+
ky
ktol
a
r
n+1
 · hn+1 .
hn+2 = min 2, max 0.6, 0.9 ·
(p)
(p−1)
ky n+1 − y n+1 k
2.2.2
Strutture in GNU Octave
Dato un generico metodo esplicito per la soluzione di un sistema differenziale (1.1), implementiamo l’algoritmo di risoluzione mediante una function [tout,yout] = method(odefun,tspan,y0), ove odefun implementa
f (t, y(t)), tspan è un vettore dei tempi ai quali si desidera la soluzione,
y0 la soluzione iniziale (corrispondente a tspan(1)), tout un vettore dei
tempi ai quali si ottiene la soluzione e yout la soluzione corrispondente.
Opzionalmente, si può passare una struttura options con i seguenti dati:
• options.RelTol tolleranza relativa
• options.AbsTol tolleranza assoluta
• options.InitialStep passo temporale iniziale
Capitolo 3
Integratori esponenziali
Il metodo Eulero esponenziale per la risoluzione di (14) è
yn+1 = exp(hA)yn + kϕ1 (hA)b(yn )
(3.1)
ove yn ≈ y(tn ), tn+1 = tn + h.
Proposizione 1. Per il metodo di Eulero esponenziale (3.1), se y(tn ) = yn ,
si ha
ky(tn+1 ) − yn+1 k = O(h2 )
e
y(tn+1 ) = yn+1 ,
se b(y(t)) = b(y0 ) ≡ b .
Dunque il metodo è esatto per problemi lineari autonomi a coefficienti costanti, di ordine 1 altrimenti.
Dimostrazione. Si ha
yn+1 = exp(hA)yn +
Z
tn+1
tn
exp((tn+1 − τ )A)g(tn )dτ ,
ove si è posto g(t) = b(y(t)). Per la formula di variazioni delle costanti (11)
Z tn+1
exp((tn+1 − τ )A)g(tn )dτ =
y(tn+1 ) − exp(hA)y(tn ) −
tn
Z tn+1
= exp(hA)y(tn ) +
exp((tn+1 − τ )A)g(τ )dτ +
tn
Z tn+1
− exp(hA)y(tn ) −
exp((tn+1 − τ )A)g(tn )dτ =
tn
Z tn+1
=
exp((tn+1 − τ )A)(g(tn ) + g ′ (τn )(τ − tn ) − g(tn ))dτ =
tn
2
= h ϕ2 (hA)g ′ (τn )
27
28
CAPITOLO 3. INTEGRATORI ESPONENZIALI
Proposizione 2. Per un problema lineare, non autonomo
(
y ′ (t) = Ay(t) + b(t), t > 0
y(t0 ) = y0
il metodo esponenziale—punto medio
yn+1 = exp(hA)yn + kϕ1 (hA)b(tn + h/2)
è di ordine 2.
Dimostrazione. Procedendo come sopra, si arriva a
Z tn+1
y(tn+1 ) = yn+1 +
exp((tn+1 − τ )A)g ′ (τn + h/2)(τ − (tn + h/2))dτ =
t
Z ntn+1
= yn+1 +
exp((tn+1 − τ )A)g ′ (τn + h/2)(τ − tn − h/2)dτ =
tn
2
= yn+1 + (h ϕ2 (hA) − h/2hϕ1 (hA))g ′ (τn + h/2) =
µ 2
¶
h I h3 A
h2 I
h3 A
4
4
= yn+1 +
+
+ O(h ) −
−
+ O(h ) g ′ (τn + h/2) .
2
6
2
2
Capitolo 4
Esercizi
1. Si consideri il seguente problema differenziale del secondo ordine ai
limiti
 ′′

 y (x) − 3 cos(y(x)) = 0, x ∈ (0, 1)
y(0) = 0

 y(1) = 1
Lo si trasformi in un sistema del primo ordine (t = x, y1 (t) = y(x),
y2 (t) = y ′ (x)) da risolvere con il metodo di Eulero esplicito e si determini, con una opportuna strategia, quale dovrebbe essere il valore
iniziale y2 (0) affinché y1 (t) = y(x) sia soluzione del problema originale.
2. Con riferimento alla Figura 4.1, l’equazione del pendolo è
 ′′

 lϑ (t) = −g sin ϑ(t)
ϑ(0) = ϑ0

 ϑ′ (0) = 0
p
La si risolva con il metodo di Crank–Nicolson fino al tempo t∗ = π l/g
(assumendo l = 1, ϑ0 = π/4). Si confronti la traiettoria con quella
del pendolo linearizzato (sin ϑ(t) ≈ ϑ(t)). Di quest’ultimo, si trovi il
numero minimo di passi temporali affinché il metodo di Eulero esplicito
produca una soluzione al tempo t∗ che dista da ϑ(t∗ ) meno di 10−2 .
3. Si calcoli y(1), ove y ′ (t) = Ay(t), y(0) = [1, . . . , 1]T , con A data da A = 100*toeplitz(sparse([-2,1,zeros(1,N-2)])), N = 10,
usando il θ-metodo con θ = 0, 1/2, 1 e diversi passi temporali h =
2−3 , 2−4 , . . . , 2−8 . Si confrontino i risultati con la soluzione di riferimento ottenuta usando θ = 1/2 e h = 2−10 , mettendo in evidenza
29
30
CAPITOLO 4. ESERCIZI
ϑ0
l
ϑ(t)
m
g
Figura 4.1: Pendolo
l’ordine del metodo usato. Si provi anche il valore θ = 2/3, discutendo
i risultati ottenuti.
4. Si risolva il sistema di ODEs
 ′

 A (t) = −2a(t)A(t)
a′ (t) = A(t)2 + Ω(t)2 − a(t)2 − 1

 Ω′ (t) = −2(a(t) + A(t))Ω(t)
con dato iniziale
(4.1)


 A(0) = 0.5
a(0) = 2

 Ω(0) = 10
con il metodo di Eulero implicito fino ad un tempo finale t∗ = 15, producendo un grafico della quantità E(t) = (A(t)2 + a(t)2 + Ω(t)2 +
1)/(2A(t)). Si confrontino le soluzioni ottenute usando 300 o 900
timesteps.
5. Si implementi il metodo Runge–Kutta–Fehlberg il cui tableaux è riportato nella Tabella 4.1, e se ne individui numericamente l’ordine. e lo si
testi sul sistema differenziale (4.1).
31
0
1
4
3
8
12
13
1
1
2
1
4
3
32
1932
2197
439
216
8
− 27
9
32
− 7200
2197
−8
7296
2197
3680
513
3544
− 2565
845
− 4104
1859
4104
− 11
40
25
216
0
1408
2565
2197
4104
− 15
16
135
0
6656
12825
28561
56430
9
− 50
2
2
55
Tabella 4.1: Metodo di Runge–Kutta–Fehlberg.
6. Si implementi la function relativa ad un generico metodo di Runge–
Kutta con tableaux dato da
c
A
bT
ove c, A e b sono dati.
Parte II
BVPs
(Problemi ai limiti)
32
33
Consideriamo il seguente problema ai limiti (boundary value problem)
 ′′
′

 u (x) = f (x, u(x), u (x)), x ∈ (a, b)
u(a) = ua
(4.2)


u(b) = ub
ove u(x) ∈ R.
Capitolo 5
Metodo di shooting
È possibile trasformare il problema (4.2) in un sistema differenziale del primo
ordine
y ′ (t) = f (t, y(t)), t ∈ (a, b]
tramite il cambiamento di variabili t = x, y1 (t) = u(x), y2 (t) = u′ (x),
f (t, y(t)) = [y2 (t), f (t, y1 (t), y2 (t)]T . Per quanto riguarda le condizioni iniziali, mentre quella per y1 (t) è y1 (a) = ua , quella per y2 (t) non è definita. Si
può allora introdurre un parametro s ∈ R e considerare la seguente famiglia
di problemi ai valori iniziali
 ′

 y (t) = f (t, y(t)), t ∈ (a, b]
y1 (a) = ua
(5.1)

 y (a) = s
2
Dato s, il sistema sopra può essere risolto con un opportuno metodo per
problemi ai valori iniziali. Poiché s è il valore della derivata prima di u(x), tale
metodo di risoluzione prende il nome di shooting. Chiamiamo y1 (t | y2 (a) = s)
la prima componente della soluzione. Si dovrà ovviamente trovare s̄ tale
che y1 (t | y2 (a) = s̄) = u(x), t = x ∈ [a, b]. In particolare, dovrà essere
y1 (b | y2 (a) = s̄) = ub . Introduciamo allora la funzione
F (s) = y1 (b | y2 (a) = s) − ub
Si tratta di risolvere l’equazione (in genere non lineare) F (s) = 0.
5.1
Metodo di bisezione
Dati due valori s1 e s2 per cui F (s1 )F (s2 ) < 0, è possibile applicare il metodo di bisezione per trovare lo zero di F (s). Poiché la soluzione di (5.1)
34
35
5.2. METODO DI NEWTON
è approssimata a meno di un errore dipendente dal passo di discretizzazione temporale, la tolleranza richiesta per il metodo di bisezione dovrà essere
(leggermente) inferiore a tale errore.
5.2
Metodo di Newton
Per applicare il metodo di Newton, è necessario calcolare F ′ (s). Definiamo
a tal scopo
v(x) =
∂
∂
u(x | u′ (a) = s) =
y1 (t | y2 (a) = s)
∂s
∂s
Derivando rispetto a s nel problema ai limiti
 ′′
′

 u (x) = f (x, u(x), u (x)),
u(a) = ua

 u′ (a) = s
x ∈ (a, b)
si ha
∂ ′′
∂
u (x) =
f (x, u(x), u′ (x))
∂s
∂s
da cui, scambiando l’ordine di derivazione
v ′′ (x) = fu (x, u(x), u′ (x))v(x) + fu′ (x, u(x), u′ (x))v ′ (x),
x ∈ (a, b)
Per quanto riguarda le condizioni iniziali per v(x), si ha
u(a | u′ (a) = s + h) − u(a | u′ (a) = s)
∂
u(a | u′ (a) = s) = lim
=
h→0
∂s
h
u a − ua
= lim
=0
h→0
h
u′ (a | u′ (a) = s + h) − u′ (a | u′ (a) = s)
∂ ′
u (a | u′ (a) = s) = lim
=
v ′ (a) =
h→0
∂s
h
s+h−s
= lim
=1
h→0
h
v(a) =
Dunque, per calcolare F ′ (s) = v(b) occorre risolvere il sistema variazionale
(lineare in v(x))
 ′′
′
′
′

 v (x) = fu (x, u(x), u (x))v(x) + fu′ (x, u(x), u (x))v (x), x ∈ (a, b)
v(a) = 0

 v ′ (a) = 1
36
CAPITOLO 5. METODO DI SHOOTING
In conclusione, per calcolare la coppia F (s) e F ′ (s) in un generico punto s,
occorre risolvere il sistema differenziale del primo ordine ai dati iniziali
 ′

 y1 (t) = y2 (t)



y2′ (t) = f (t, y1 (t), y2 (t))





y3′ (t) = y4 (t)



 y ′ (t) = f (t, y (t), y (t))y (t) + f (t, y (t), y (t))y (t)
y1
1
2
3
y2
1
2
4
4

y1 (a) = ua





y2 (a) = s




y3 (a) = 0




y4 (a) = 1
fino al tempo t = b. Quindi F (s) = y1 (b) e F ′ (s) = y3 (b). Poiché le equazioni
per y1′ (t) e y2′ (t) non dipendono da y3 (t) e y4 (t), è possibile disaccoppiare le
prime due componenti dalle seconde due.
Una semplificazione del metodo di Newton che non richiede il calcolo di
F ′ (s) è il metodo delle secanti.
5.3
Problema ai limiti con frontiera libera
Un caso particolarmente interessante per l’applicazione del metodo di shooting è quello a frontiera libera (free boundary)
 ′′
u (x) = f (x, u(x), u′ (x)),




u(s) = α

u′ (s) = β



u(b) = ub
x ∈ (s, b)
(5.2)
ove i valori di u e di u′ sono assegnati in un punto incognito s, s < b. La
funzione di cui si deve trovare lo zero è, in questo caso,
F (s) = u(b | u(s) = α, u′ (s) = β) − ub
(scriveremo F (s) = u(b | s) − ub per brevità). Dati due punti s1 e s2 tali
che F (s1 )F (s2 ) < 0, l’applicazione del metodo di bisezione non presenta
difficoltà. Per quanto riguarda il metodo di Newton, il sistema variazionale
per
∂
u(x | s + h) − u(x | s)
v(x) =
u(x | s) = lim
h→0
∂s
h
5.3. PROBLEMA AI LIMITI CON FRONTIERA LIBERA
37
è analogo al caso precedente. L’unica diversità è data dalle condizioni iniziali
(in s). Si ha
u(s | s + h) − u(s | s)
v(s) = lim
h→0
h
Ora, u(s | s) = α. Poi
u(s | s + h) = u(s + h | s + h) − hu′ (s + h | s + h) + O(h2 ) = α − hβ + O(h2 )
Dunque, v(s) = −β. In maniera analoga
u′ (s | s + h) − u′ (s | s)
= −u′′ (s)
h→0
h
v ′ (s) = lim
ove il valore u′′ (s) si ricava dal problema (5.2).
Capitolo 6
Differenze finite
Consideriamo il seguente problema ai limiti
u′′ (x) = f (x, u(x), u′ (x)),
x ∈ (a, b)
(6.1)
soggetto a condizioni al bordo di Dirichlet o di Neumann.
6.1
Costruzione delle matrici delle differenze
finite
La matrice corrispondente all’approssimazione mediante differenze finite di
ordine 2 della derivata prima con griglia non equispaziata è


 

′
∗
∗
∗
∗
∗
∗
u (x1 )
u(x
)
1
1
−1

0
0
...
0
 u′ (x2 )  
  u(x2 ) 
h1 +h2
h1 +h2

 



.
...
 u′ (x3 )   0

1
−1
.

u(x
)
0
.
3

 


h2 +h3
h2 +h3

≈

  .


..
..

.
.
.
.




.
.
.
.
.
.
.
.
.
.

0
.
 ′ .
 



−1
1
u (xm−1 )  0
u(xm−1 )

...
0
0
hm−2 +hm−1
hm−2 +hm−1
u′ (xm )
u(xm )
∗
∗
∗
∗
∗
∗
Dati i nodi x (vettore colonna di lunghezza m), è possibile costruire il vettore
[h1 , h2 , . . . , hm−1 ]T con il comando h=diff(x). Allora la matrice, a meno
della prima e dell’ultima riga, può essere costruita, direttamente in formato
sparso, con i comandi
> d = 1./(h(1:m-2)+h(2:m-1));
> spdiags([[-d;0;0],[0;0;d]],[-1,1],m,m)
38
6.2. CONDIZIONI DI DIRICHLET
39
La costruzione della matrice relativa alla derivata seconda è analoga. Nel caso
di griglia equispaziata, di passo h, le matrici relative alle approssimazione
della derivata prima e seconda possono essere costruite con i comandi
> toeplitz(sparse([0,-1/(2*h),zeros(1,m-2)]),...
> sparse([0,1/(2*h),zeros(1,m-2)]))
e
> toeplitz(sparse([-2/h^2,1/h^2,zeros(1,m-2)]))
rispettivamente.
6.2
Condizioni di Dirichlet
Se vengono prescritti i valori u(a) = ua o u(b) = ub , conviene discretizzare,
in un primo momento, il problema ai limiti senza tener conto delle condizioni
al bordo. Per esempio, la discretizzazione del problema ai limiti
 ′′

 u (x) = 0, x ∈ (a, b)
u(a) = ua

 u(b) = u
b
diventa

−2
1

1 
0

h2  ...

0
0
  

0
u1
1
0 0 ··· 0




−2 1 0 · · · 0   u2  0

. . . . . . . . . . . . ..   ..   .. 
.   .  .
  .  = .
... ... ... ...
 .  .
0
  .  .
· · · 0 1 −2 1  um−1  0
0
··· 0 0
1 −2
um
Poi, si correggono le equazioni relative ai nodi al bordo

1
1

1 
0

2
h  ...

0
0
 


0
0 0 ··· 0
ua
u1





−2 1 0 · · · 0  u2 
 0. 
.
.
... ... ... ... .  . 

1 
.  . 
 .. 
=




. . . . . . . . . . . .   ..  h2  .. 

0  . 
.
0
· · · 0 1 −2 1 um−1 
··· 0 0
0 1
ub
um
40
6.3
CAPITOLO 6. DIFFERENZE FINITE
Condizioni di Neumann
Se vengono prescritti i valori della derivata prima u′ (a) = u′a o u′ (b) = u′b ,
è necessario approssimare la derivata prima con uno stencil non simmetrico, eventualmente con ordine di approssimazione minore. Per esempio, la
discretizzazione del problema ai limiti
 ′′

 u (x) = 0, x ∈ (a, b)
u′ (a) = u′a

 u(b) = u
b
potrebbe essere

−1
1

1 
0

2
h  ...

0
0
o, meglio,

  ′ 
ua /h
1
0 0 ··· 0
u1
 u2   0 
−2 1 0 · · · 0

 

. . . . . . . . . . . . ..   u   .. 
.  3   . 
 ..  =  . 
... ... ... ... 

  . 
0
 .   . 
· · · 0 1 −2 1 um−1   0 
um
ub /h2
··· 0 0
0 1

  ′ 
 3
ua /h
− 2 2 − 21 0 · · · 0
u1




 1 −2 1
0 · · · 0  u 2   0 


. . . . . . . . . . . . ..   u   .. 
1 
3
. 
  . 
 0
 ..  =  . 

... ... ... ... 


 . 
h2  ...
0  . 
  . 

 0 ··· 0
1 −2 1 um−1   0 
um
0 ··· 0
0
0 1
ub /h2
Infatti lo stencil in avanti [−3, 4, −1]/(2h) produce una approssimazione del
secondo ordine in h della derivata prima.
6.4
Differenze finite centrate del secondo ordine
Sia u ∈ C 4 ([a, b]) e xi = a + (i − 1)h, i = 1, 2, . . . , m, h = (b − a)/(m − 1).
Sviluppando in serie, si ha
h2 ′′
h3
h4
u (xi ) + u(3) (xi ) + u(4) (x̂i )
2
6
24
3
2
h
h4
h
u(xi−1 ) = u(xi ) − hu′ (xi ) + u′′ (xi ) − u(3) (xi ) + u(4) (x̃i )
2
6
24
u(xi+1 ) = u(xi ) + hu′ (xi ) +
41
6.5. CONVERGENZA PER UN PROBLEMA MODELLO
da cui
u(xi+1 ) − 2u(xi ) + u(xi−1 )
+ τi
h2
2
ove τi = − h12 u(4) (x̄i ) è l’errore locale (u(4) (x̂i )/24 + u(4) (x̃i )/24 = u(4) (x̄i )/12,
per un opportuno x̄i , per il teorema dei valori intermedi).
u′′ (xi ) =
6.5
Convergenza per un problema modello
Consideriamo il seguente problema modello (elasticità della trave)

′′

 −u (x) + q(x)u(x) = g(x), x ∈ (a, b)
u(a) = ua

 u(b) = u
b
(6.2)
con q, g ∈ C 0 ([a, b]), q(x) ≥ 0 per x ∈ [a, b]. Vogliamo studiare l’esistenza,
l’unicità e la regolarità della soluzione analitica.
6.5.1
Unicità
Siano u1 (x) e u2 (x) due soluzioni di (6.2). Allora, z(x) = u1 (x) − u2 (x)
soddisfa il problema omogeneo
 ′′

 −z (x) + q(x)z(x) = 0, x ∈ (a, b)
z(a) = 0
(6.3)

 z(b) = 0
Proposizione 3. Se z(x) è soluzione di (6.3), allora z(x) ≡ 0.
Dimostrazione (metodo dell’energia). Moltiplicando l’equazione per z(x) ed
integrando si ha
Z b
Z b
′′
q(x)z(x)2 dx =
−z (x)z(x)dx +
0=
a
a
Z b
Z b
′
b
′
2
= [−z (x)z(x)]a +
z (x) dx +
q(x)z(x)2 dx =
a
a
Z b
Z b
q(x)z(x)2 dx
z ′ (x)2 dx +
=
a
a
Poiché le funzioni integrande sono non negative, si ha che deve essere necessariamente z(x) ≡ 0.
Dunque, u1 (x) ≡ u2 (x).
42
6.5.2
CAPITOLO 6. DIFFERENZE FINITE
Esistenza
Sia z(x) = c1 z1 (x) + c2 z2 (x) la soluzione generale di −z ′′ (x) + q(x)z(x) = 0.
La soluzione di (6.3) si ottiene imponendo
(
c1 z1 (a) + c2 z2 (a) = 0
c1 z1 (b) + c2 z2 (b) = 0
Poiché sappiamo che z(x) ≡ 0 è l’unica soluzione, si ha che la matrice
¸
·
z1 (a) z2 (a)
z1 (b) z2 (b)
è non singolare.
La soluzione generale di −u′′ (x) + q(x)u(x) = g(x) è u(x) = c1 z1 (x) +
c2 z2 (x) + s(x) (s(x) soluzione particolare). La soluzione di (6.2) si ottiene
imponendo
(
c1 z1 (a) + c2 z2 (a) = ua − s(a)
c1 z1 (b) + c2 z2 (b) = ub − s(b)
cioè risolvendo un sistema lineare non singolare che ammette dunque (unica)
soluzione.
Si è costretti a ridursi ad un problema modello perché problemi ai limiti
anche molto semplici possono non avere soluzione: si consideri, per esempio,
 ′′

 u (x) + u(x) = 0
u(0) = 0

 u(π) = 1
La soluzione generale è c1 cos(x) + c2 sin(x), ma non è possibile imporre le
condizioni al bordo.
6.5.3
Regolarità
Proposizione 4. Se q, g ∈ C k ([a, b]), allora u ∈ C k+2 ([a, b]).
Dimostrazione. Se q, g ∈ C 0 ([a, b]), poiché la soluzione u esiste, u′′ è definita
in ogni punto x ∈ [a, b], e dunque u′ esiste (ed è derivabile). Quindi u ∈
C 0 ([a, b]) e quindi u′′ ∈ C 0 ([a, b]). Dunque u ∈ C 2 ([a, b]). Sia vero adesso
l’enunciato per k e siano q, g ∈ C k+1 ([a, b]): poiché anche u ∈ C k+1 ([a, b]), si
ha u′′ ∈ C k+1 ([a, b]) da cui u ∈ C k+3 ([a, b]).
6.5. CONVERGENZA PER UN PROBLEMA MODELLO
43
Ci occupiamo adesso di analizzare la convergenza del problema discretizzato mediante differenze finite centrate del secondo ordine, che si scrive

ui+1 − 2ui + ui−1


+ qi ui = gi , i = 2, 3, . . . , m − 1
−
h2
u1 = u a



um = u b
ove qi = q(xi ) e gi = g(xi ).
6.5.4
Consistenza
Se si sostituisce ui con la soluzione analitica u(xi ), si ottiene

u(xi+1 ) − 2u(xi ) + u(xi−1 )


+ q(xi )u(xi ) − g(xi ) = τi ,
−
h2
i = 2, 3, . . . , m − 1
u(x1 ) = ua



u(xm ) = ub
da cui si deduce che il metodo numerico è consistente di ordine 2. Il sistema
lineare da risolvere per trovare uh = [u1 , u2 , . . . , um−1 , um ]T è
1
0
...
−1 2 + q2 h2
−1


1  0
−1
2 + q3 h2

...
...
h2  ...

 0
...
0
0
...
...
...
0
...
...
−1
...
0
...
−1
0
2 + qm−1 h2
0
0
...
...
−1
..
.
..
.
0
..
.
..
.
...
..
.
..
.

 


0
u1
ua /h2
 u2   g2 
0

 

.. 
 


.   u3   g3 
  ..  =  .. 
 


0
 .   . 
−1 um−1   gm−1 
um
ub /h2
1
e può essere semplificato in

2 + q2 h2
−1

 −1
2 + q3 h2


..
.
1 
0

.
h2 
..

0


..

.
...
0
...
0
...
−1
0
2 + qm−2 h2
−1
cioè
Ah uh = gh
0
..
.


 

g2 + ua /h2
 u2
  u3  

g3
 


  ..  

..





.
0
 .  = 

  ..  

..





.
.
0
 


 um−2  

gm−2

−1
2
u
g
+
u
/h
m−1
m−1
b
2 + qm−1 h2
(6.4)
44
6.5.5
CAPITOLO 6. DIFFERENZE FINITE
Esistenza ed unicità
Proposizione 5. Il sistema lineare (6.4) è non singolare e dunque ammette
un’unica soluzione.
Dimostrazione (metodo dell’energia discreto). Consideriamo il prodotto z T Ah z
ove z = [z2 , z3 , . . . , zm−1 ]T . Si ha
1
[(2 + q2 h2 )z22 − z2 z3 − z3 z2 + (2 + q3 h2 )z32 − z3 z4 + . . . +
2
h
2
+ . . . − zm−1 zm−2 + (2 + qm−1 h2 )zm−1
]=
1
2
]+
= 2 [z22 + (z2 − z3 )2 + (z3 − z4 )2 + . . . + (zm−2 − zm−1 )2 + zm−1
h
m−1
X
+
qi zi2 ≥ 0
z T Ah z =
i=2
Dunque la matrice Ah è definita positiva e quindi non singolare.
6.5.6
Proprietà di Ah
Ah è una matrice simmetrica e diagonalmente dominante. È possibile usare
i metodi iterativi, semi-iterativi e diretti senza pivoting per la soluzione del
sistema lineare. Inoltre, è una M -matrice, cioè i suoi elementi extra-diagonali
sono non positivi e la sua inversa ha elementi non negativi.
6.5.7
Stabilità
Consideriamo due soluzioni relative a dati perturbati g̃h e ĝh . Si ha
Ah ũh = g̃h
Ah ûh = ĝh
da cui
(ũh − ûh ) = A−1
h (g̃h − ĝh )
Se si vuole che le perturbazioni sui dati non si ripercuotano in maniera distruttiva sulle soluzioni, occorre che la matrice A−1
h sia limitata in norma
indipendentemente da h, in particolare per h → 0. Consideriamo la matrice
Ah,q=0 corrispondente alla stessa discretizzazione nel caso q(x) ≡ 0. Si ha
Ah − Ah,q=0 = diag(q2 , . . . , qm−1 ) ≥ 0. Allora
−1
−1
−1
A−1
h,q=0 − Ah = Ah,q=0 (Ah − Ah,q=0 )Ah ≥ 0
6.5. CONVERGENZA PER UN PROBLEMA MODELLO
45
−1
perché Ahq =0 e Ah sono M -matrici. Allora A−1
h ≤ Ah,q=0 . Osserviamo poi
T
che A−1
0 [1, . . . , 1] è la soluzione discreta (approssimata) di
 ′′

 −v (x) = 1
v(a) = 0

 v(b) = 0
la cui soluzione analitica è v(x) = (x − a)(b − x)/2. Poiché v (3) (x) ≡ 0, cosı̀
è per v (4) e dunque l’errore locale, per questo problema, è 0. Dunque
−1
T
kA−1
h,q=0 k∞ = kAh,q=0 [1, . . . , 1] k∞ =
=
max
i=2,...,m−1
max
i=2,...,m−1
v(xi ) ≤ max v(x) ≤
x∈[a,b]
vi =
(b − a)2
8
e poiché kA−1
h k∞ ≤ kAh,q=0 k∞ , si ha la maggiorazione richiesta.
6.5.8
Convergenza
Definiamo eh = [e2,h , . . . , em−1,h ]T = [u2,h − u(x2 ), . . . , um−1,h − u(xm−1 )]T .
Poiché
Ah [u2,h , . . . , um−1,h ]T = gh
Ah [u(x2 ), . . . , u(xm−1 )]T = gh + τh
si deduce eh = −A−1
h τh . Combinando i risultati di consistenza e stabilità, si
ottiene, per il problema (6.2) discretizzato mediante differenze finite centrate
del secondo ordine,
(b − a)2 h2 (4)
ku k∞
keh k∞ ≤
8
12
e dunque l’errore è proporzionale a h2 , posto che u ∈ C 4 ([a, b]).
6.5.9
Un esempio: l’equazione della catenaria
Consideriamo l’equazione della catenaria

p
′′
′
2

 u (x) = a 1 + u (x) ,
u(−1) = 1


u(1) = 1
x ∈ (−1, 1)
46
CAPITOLO 6. DIFFERENZE FINITE
La discretizzazione mediante differenze finite centrate del secondo ordine è




u1
1
¡ u3 −u1 ¢2 
 u2 
 q

 . 

1
+
2h

 . 

.




..
A .  − a
=b
.

 .. 
q
¢
¡



um −um−2 2 

um−1 
 1+
2h
1
um
Si tratta dunque di risolvere il sistema non lineare
p
F (u) = Au − a 1 + (Bu)2 − b = 0
ove

1 0
0
1 −2 1


.
.
1 
0 .. ..
A= 2
.
h 
 .. . . . . . .

0 · · · 0
0 ··· 0
0 ···
0 ···
.. ..
.
.
.. ..
.
.
1 −2
0
0
Lo jacobiano di F (u) è
J(u) = A − aD(u)B,
6.6
Norme

0
0

.. 
.
,

0

1
1

0
0
−1 0


..
.
1 
0
B=

.
.
2h  ..
..

 0 ···
0 ···
D = (dij (u)),
0
0 ···
1
0 ···
.. ..
..
.
.
.
.. ..
..
.
.
.
0 −1 0
0
0
0

(Bu)i

q
,
1 + (Bu)2i
dij (u) =


0,

0
0

.. 
.
,

0

1
0

1
h2

−a
 0 




b =  ... 


 0 
1
−a
h2
i=j
i 6= j
Data una funzione u(x) e due discretizzazioni su nodi equispaziati [ũ1 , ũ2 , . . . , ũm ] ≈
[u(x̃1 ), u(x̃2 ), . . . , u(x̃m )] e [û1 , û2 , . . . , ûl ] ≈ [u(x̂1 ), u(x̂2 ), . . . , u(x̂l )], {x̃i }i ⊂ [a, b],
{x̂i }i ⊂ [a, b], con m 6= l, non ha molto senso confrontare gli errori k[u(x̃1 ) −
ũ1 , u(x̃2 )−ũ2 , . . . , u(x̃m )−ũm ]k2 e k[u(x̂1 )−û1 , y(x̂2 )−û2 , . . . , u(x̂l )−ûl ]k2q
. Si preferisce usare la norma infinito, oppure la norma k[v1 , v2 , . . . , vm ]kL2 = kvk2 b−a
m , che
risulta essere una approssimazione mediante quadratura con formula dei rettangoli
della norma in L2 .
Capitolo 7
Equazione di Poisson
Di particolare interesse è l’equazione di Poisson
−∇2 u(x) = f (x),
x ∈ Ω ⊂ Rd
ove ∇2 è l’operatore laplaciano definito da
∇2 =
d
X
∂2
∂x2k
k=1
L’equazione è solitamente accompagnata da condizioni al bordo di Dirichlet o di
Neumann.
7.1
Equazione di Poisson bidimensionale
Analizziamo numericamente in dettaglio il caso d = 2 (x = (x, y)) e Ω = [a, b] ×
[c, d].
7.1.1
Condizioni al bordo di Dirichlet
Consideriamo dapprima il caso con condizioni al bordo di Dirichlet. Dunque

−∇2 u(x, y) = f (x, y), (x, y) ∈ [a, b] × [c, d] ⊂ R2






 u(a, y) = Da (y)
u(b, y) = Db (y)



u(x, c) = Dc (x)




u(x, d) = Dd (x)
con le necessarie condizioni di compatibilità ai vertici. Introduciamo una discretizzazione xj = a+(j−1)hx , j = 1, 2, . . . , mx , hx = (b−a)/(mx −1) e yi = c+(i−1)hy ,
47
48
CAPITOLO 7. EQUAZIONE DI POISSON
i = 1, 2, . . . , my , hy = (d − c)/(my − 1). Introduciamo infine la discretizzazione di
u(x, y) definita da
uk ≈ u(xj , yi ), k = (j − 1)my + i
di cui si vede un esempio in Figura 7.1. La matrice di discretizzazione alle dif-
u4
u8
u12
u16
u3
u7
u11
u15
u2
u6
u10
u14
u1
u5
u9
u13
Figura 7.1: Numerazione di una griglia bidimensionale
ferenze finite centrate del secondo ordine, senza tener conto delle condizioni al
bordo, è data da
A = Imx ⊗ Ay + Ax ⊗ Imy
ove ⊗ indica il prodotto di Kronecker e
2
−1
0
6
6−1
6
6
6
1 6
0
Ay = 2 6
hy 6
6
6 0
6
6 .
4 ..
0
2
..
.
−1
..
.
..
2
..
...
...
...
...
.
0
..
.
..
...
..
.
..
.
0
...
−1
0
.
.
2
−1
3
0
.7
.7
.7
7
7
07
7,
7
7
07
7
7
−15
2
2
2
6
6−1
6
6
6
1 6
0
Ax = 2 6
hx 6
6
6 0
6
6 .
4 ..
0
0
...
...
2
..
.
−1
..
.
0
..
.
..
..
.
..
...
..
.
..
.
0
...
−1
0
−1
.
...
...
.
2
−1
3
0
.7
.7
.7
7
7
07
7
7
7
07
7
7
−15
2
ove Ay ∈ Rmy ×my e Ax ∈ Rmx ×mx . Poi, le righe di indice, diciamo k, corrispondente ad un nodo al bordo vanno sostituite con il vettore della base canonica
ek , diviso per h2x + h2y . Il termine noto è [b1 , b2 , . . . , bmx my ]T , ove


f (xj , yi ) se (xj , yi ) è un nodo interno, k = (j − 1)my + i



Da (yi )


se xj = a, k = (j − 1)my + i

2
2

y
 hDx +h
(y
)
b i
se xj = b, k = (j − 1)my + i
bk = h2x +h2y


D
(x
)
c
j

se yi = c, k = (j − 1)my + i
 h2 +h2

x
y


D
(x
)

j
d
 2 2
se yi = d, k = (j − 1)my + i
h +h
x
y
49
7.1. EQUAZIONE DI POISSON BIDIMENSIONALE
Alternativamente, si può sostituire il solo termine diagonale delle righe corrispondenti ad un nodo al bordo con un coefficiente M/(h2x + h2y ), M ≫ 1 e moltiplicare
per M il corrispondente elemento nel termine noto. Questa procedura permette
di assegnare, di fatto, le condizioni al bordo di Dirichlet, mantenendo la matrice
A simmetrica.
7.1.2
Condizioni al bordo miste
L’equazione di Poisson non può essere accompagnata solo da condizioni al bordo di
Neumann, altrimenti la soluzione è indeterminata. Consideriamo allora il seguente
problema con condizioni al bordo miste

−∇2 u(x, y) = f (x, y), (x, y) ∈ [a, b] × [c, d] ⊂ R2






 u(b, y) = Db (y)


 u(x, c) = Dc (x), Dc (b) = Db (c)
∂u


− (x, y) = Na (y), x = a, c < y < d


∂x




∂u


(x, y) = Nd (x), y = d, x < b
∂y
La matrice di discretizzazione alle differenze finite centrate del secondo ordine è
data da
0
0
⊗ Ay + Ax ⊗ Im
A = Im
y
x
ove
2
−1
0
6
6−1
6
6
6
1 6
0
Ay = 2 6
hy 6
6
6 0
6
6 .
4 ..
0
2
..
.
−1
..
.
..
2
e
0
Im
x
..
...
...
...
...
.
0
..
.
..
...
..
.
..
.
0
...
−1
1
2
.
.
2
−2
3
0
.7
.7
.7
7
7
07
7,
7
7
07
7
7
−1 5
3
2


0 0
0 ... 0

.
0 1
0 . . . .. 



 .. .. ..
,
= 0
.
.
. 0



 ..
. . . . 0
1 0
0 ... ... 0 1
2
3
2
6
6 −1
6
6
6
1 6
0
Ax = 2 6
hx 6
6
6 0
6
6 .
4 ..
0
0
Im
y
1
2
...
...
2
..
.
−1
..
.
0
..
.
..
..
.
..
...
..
.
..
.
0
...
−1
0
−2
.
...
...


1 0
0 ... 0

.
0 1
0 . . . .. 



 .. .. ..
= 0
.
.
. 0



 ..
. . . . 0
1 0
0 ... ... 0 0
.
2
−1
3
0
.7
.7
.7
7
7
07
7
7
7
07
7
7
−15
2
Poi, le righe di indice, diciamo k, corrispondente ad un nodo al bordo su cui
sono prescritte condizioni di Dirichlet vanno sostituite con il vettore della base
canonica ek , diviso per h2x + h2y . La riga di indice my , corrispondente al nodo di
50
CAPITOLO 7. EQUAZIONE DI POISSON
bordo (a, c), va sostituita con
[1, 0, . . . , 0] ⊗
1
[0, . . . , 0, 1/2, −2, 3/2]
h2y
Il termine noto è [b1 , b2 , . . . , bmx my ]T , ove
bk =


f (xj , yi )



Db (yi )



 h2x +h2y
Dc (xj )
h2x +h2y


N

a (yi )


hx


 Nd (xj )
hy
se (xj , yi ) è un nodo interno, k = (j − 1)my + i
se xj = b, k = (j − 1)my + i
se yi = c, k = (j − 1)my + i
se xj = a, k = (j − 1)my + i, i 6= 1, i 6= my
se yi = d, k = (j − 1)my + i, j 6= mx
Capitolo 8
Metodi variazionali
8.1
Formulazione variazionale di un problema modello
Consideriamo il seguente problema ai limiti (equazione di Poisson)
(
−u′′ (x) = f (x), x ∈ (0, 1)
u(0) = u(1) = 0
(8.1)
con f ∈ C 0 ([0, 1]). Introduciamo il seguente spazio lineare:
V = {v: v ∈ C 0 ([0, 1]), v ′ continua a tratti e limitata, v(0) = v(1) = 0}
e il funzionale lineare J : V → R dato da
1
J(v) = (v ′ , v ′ ) − (f, v)
2
ove
(v, w) =
Z
1
v(x)w(x)dx
0
Teorema 2 (Formulazione variazionale). Se u è la soluzione del problema (8.1),
allora
(u′ , v ′ ) = (f, v), ∀v ∈ V
(8.2)
e, equivalentemente,
J(u) ≤ J(v),
∀v ∈ V
Dimostrazione. Sia u soluzione di (8.1). Allora, per ogni v ∈ V ,
Z
0
1
′′
−u (x)v(x)dx =
Z
1
f (x)v(x)dx = (f, v)
0
51
(8.3)
52
CAPITOLO 8. METODI VARIAZIONALI
Integrando per parti,
Z 1
¯1 Z
¯
′′
′
−u (x)v(x)dx = −u (x)v(x)¯ +
0
0
1
u′ (x)v ′ (x)dx = (u′ , v ′ )
0
poiché v(0) = v(1) = 0.
Sia adesso u ∈ V soluzione di (8.2) e w = v − u, per v ∈ V . Allora w ∈ V e
1
J(v) = J(u + w) = (u′ + w′ , u′ + w′ ) − (f, u + w) =
2
1 ′ ′
1
= (u , u ) + (u′ , w′ ) + (w′ , w′ ) − (f, u) − (f, w) ≥ J(u)
2
2
′
′
′
′
perché (u , w ) − (f, w) = 0 e (w , w ) ≥ 0. Dunque J(u) ≤ J(v).
Sia infine u ∈ V soluzione di (8.3). Allora
J(u) ≤ J(u + εv),
∀ε, ∀v ∈ V
Allora ψ(ε) = J(u + εv) ha un minimo in ε = 0 e dunque ψ ′ (0) = 0. Poiché
ε2
1
ψ(ε) = (u′ , u′ ) + ε(u′ , v ′ ) + (v ′ , v ′ ) − (f, u) − ε(f, v)
2
2
si ha
i
h
ε
ψ(ε) − ψ(0)
= lim (u′ , v ′ ) + (v ′ , v ′ ) − (f, v) =
ε→0
ε→0
ε
2
= (u′ , v ′ ) − (f, v)
0 =ψ ′ (0) = lim
Abbiamo dunque dimostrato le seguenti implicazioni
(8.1) ⇒ (8.2) ⇔ (8.3)
Verifichiamo che la soluzione di (8.2) è unica: se u1 e u2 sono due soluzioni, allora
(u′1 − u′2 , v ′ ) = 0,
∀v ∈ V
e in particolare per v = u1 − u2 . Dunque
Z 1
(u′1 (x) − u′2 (x))2 dx = 0
0
u′1 (x)
u′2 (x)
e quindi
−
= (u1 (x) − u2 (x))′ = 0. Pertanto u1 − u2 è costante e
siccome u1 (0) − u2 (0) = 0, allora u1 (x) − u2 (x) = 0.
La soluzione di (8.1) si chiama soluzione forte del problema (8.1), mentre la
soluzione di (8.2) (o, equivalentemente, di (8.3)) si chiama soluzione debole del
problema (8.1). Perché soluzione debole? Se u è soluzione di (8.2) e u ∈ C 2 ([0, 1]),
allora 0 = (u′ − f, v) = (−u′′ − f, v) per ogni v ∈ V . Poiché u′′ + f è continua, si
deduce −u′′ (x) = f (x) per 0 < x < 1.
Le formulazioni variazionali del problema (8.1) sono in realtà le più “fisiche”:
pensando al problema della trave, esse permettono di descrivere anche il caso in
cui il carico f (x) non sia continuo (ma, per esempio, applicato in un solo punto).
8.1. FORMULAZIONE VARIAZIONALE DI UN PROBLEMA MODELLO53
8.1.1
Metodo di approssimazione variazionale
Prendiamo un sottospazio Vh di V di dimensione finita. Si cerca allora uh ∈ Vh
tale che
(uh , vh )′ = (f, vh ), ∀vh ∈ Vh
(8.4)
(metodo di Galerkin) ove, per brevità, (uh , vh )′ = (u′h , vh′ ), o, equivalentemente
J(uh ) = inf J(vh )
vh ∈Vh
(metodo di Ritz).
Teorema 3. Il problema (8.4) ha un’unica soluzione.
Dimostrazione. Sia {φj }m
j=1 una base di Vh . Allora
uh (x) =
m
X
uj φj (x)
j=1
e il problema (8.4) si riscrive, per i = 1, 2, . . . , m,
Z
1
0

′
m
m
X
X
u′h (x)φ′i (x)dx = 
uj φj , φi  =
(φj , φi )′ uj = Auh = (f, φi )
j=1
j=1
ove A = (φj , φj )′ . Se (uh , vh )′ = 0 per ogni vh ∈ Vh (cioè uh è soluzione del
problema omogenero), allora, in particolare, (uh , uh )′ = (u′h , u′h ) = 0 da cui uh = 0.
D’altra parte, se (uh , vh )′ = 0 per ogni vh ∈ Vh , allora uh = 0 è soluzione di
m
X
(φj , φi )′ uj = A[u1 , u2 , . . . , um ]T = 0,
i = 1, 2, . . . , m
j=1
Dunque, la matrice A = (φj , φi )′ = (φ′j , φ′i ) è non singolare.
La matrice A, che risulta essere simmetrica e definita positiva, si chiama matrice di rigidezza (stiffness matrix) e il vettore (f, φi ) vettore di carico (load vector).
Vale poi il seguente risultato:
Teorema 4. Se u è soluzione di (8.2) e uh di (8.4), allora
ku − uh k′ ≤ inf ku − vh k′
vh ∈Vh
ove k·k′ =
p
(·, ·)′ .
54
CAPITOLO 8. METODI VARIAZIONALI
Dimostrazione. Dalle uguaglianze
(u, v)′ = (f, v)
∀v ∈ V e, dunque, ∀v ∈ Vh
′
(uh , vh ) = (f, vh ) ∀vh ∈ Vh
si ricava ((u − uh ), wh )′ = 0 per ogni wh ∈ Vh . Dunque
(u − uh , u − uh )′ = (u − uh , u − vh + vh − uh )′ = (u − uh , u − vh )′ ≤
≤ ku − uh k′ ku − vh k′
(per la disuguaglianza di Cauchy–Schwartz) da cui
ku − uh k′ ≤ ku − vh k′ ,
∀vh ∈ Vh
e quindi la tesi.
Per definizione, uh è allora la proiezione ortogonale della soluzione esatta u sul
sottospazio Vh , tramite il prodotto scalare (·, ·)′ .
La scelta di Vh caratterizza il metodo. Da un lato bisogna considerare la
regolarità della soluzione richiesta. Dall’altro la difficoltà di assemblare la matrice
di rigidezza e di risolvere il sistema lineare.
Stabilità e consistenza
La consistenza del metodo di Galerkin discende da
(u, vh )′ = (f, vh ),
∀vh ∈ Vh
(il metodo si dice fortemente consistente). Per quanto riguarda la stabilità, se
dimostriamo che
p
(8.5)
kuh k′ ≤ C (f, f )
allora abbiamo
q
kũh − ûh k ≤ C (f˜ − fˆ, f˜ − fˆ)
′
e cioè che piccole variazioni sui dati producono piccole variazioni sulle soluzioni.
Per dimostrare (8.5), osserviamo che
Z 1
Z 1
¯1 Z 1
′
2
2 ¯
′
uh (x)2 dx
xuh (x) dx = uh (x) x¯ −
2xuh (x)uh (x)dx =
da cui
e quindi
0
0
0
0
q
p
p
p
(uh , uh ) ≤ 2 (uh , uh ) (u′h , u′h ) = 2 (uh , uh ) (uh , uh )′
Poiché uh soddisfa
p
(uh , uh ) ≤ 2kuh k′
(uh , uh )′ = (f, uh )
si ricava, suppondendo f a quadrato sommabile,
p
p
p
kuh k′2 ≤ (f, f ) (uh , uh ) ≤ 2 (f, f )kuh k′
da cui la tesi.
8.1. FORMULAZIONE VARIAZIONALE DI UN PROBLEMA MODELLO55
Metodo degli elementi finiti (FEM)
Introduciamo una discretizzazione dell’intervallo [0, 1] a passo variabile, come in
m−1
Figura 8.1 Lo spazio Vh è generato dalle funzioni di base {φj }j=2
, le quali sono
hj−1
h1
x1
φj
φj−1
φ1
x2
xj−2
xj−1
xj
φj+1
hm−1
hj
xj+1
xm
Figura 8.1: Hat functions
definite da
e
mentre, per j = 1,
e
e, per j = m,
e

x − xj−1


, xj−1 ≤ x ≤ xj



 hj−1
φj (x) = xj+1 − x , x ≤ x ≤ x
j
j+1


hj



 0,
altrimenti

1


, xj−1 < x < xj



 hj−1
1
φ′j (x) =
− , xj < x < xj+1


hj




0,
altrimenti

 x2 − x , x1 ≤ x ≤ x2
h1
φ1 (x) =

0,
altrimenti

 − 1 , x1 < x < x2
′
h1
φ1 (x) =

0,
altrimenti
x−x
m−1

, xm−1 ≤ x ≤ xm
h
m−1
φm (x) =

0,
altrimenti
φ′m (x)
=


1
hm−1

0,
, xm−1 < x < xm
altrimenti
56
CAPITOLO 8. METODI VARIAZIONALI
Dunque, nell’approssimazione
uh (x) =
m
X
uj φj (x)
j=1
i coefficienti uj sono i valori di uh nei nodi xj . Per 1 < j < m,
µ
¶2
¶
Z xj
Z xj +hj µ
1
1
1
1 2
′
ajj = (φj , φj ) =
+
dx +
dx =
−
hj−1
hj
hj−1 hj
xj −hj−1
xj
Z xj +hj
1
1 1
− · dx = −
aj j+1 = aj+1 j = (φj , φj+1 )′ =
hj hj
hj
xj
Per j = 1 e j = m, si ha invece
¶
Z x1 +h1 µ
1
1 2
dx =
−
a11 =
h1
h1
x1
Z x1 +h1
1 1
1
a12 = a21 =
−
dx = −
h
h
h
1 1
1
x1
Z xm
1
1
1
dx = −
am−1 m = am m−1 =
−
hm−1 hm−1
hm−1
xm −hm−1
µ
¶2
Z xm
1
1
amm =
−
dx =
hm−1
hm−1
xm −hm−1
Per quanto riguarda il calcolo di (f, φi ) si può ricorrere alla formula di quadratura
del trapezio che risulta essere sufficientemente accurata. Si ha dunque, per 1 <
i < m,
Z xi
Z xi+1
x − xi−1
xi+1 − x
fi = (f, φi ) =
f (x)
dx +
dx ≈
f (x)
hi−1
hi
xi−1
xi
hi−1
hi
hi−1 + hi
≈ f (xi )
+ f (xi ) = f (xi )
2
2
2
Per i = 1 e i = m si ha invece
Z x2
x2 − x
h1
f1 = (f, φ1 ) =
f (x)
dx ≈ f (x1 )
h1
2
Zx1xm
hm−1
x − xm−1
dx ≈ f (xm )
fm = (f, φm ) =
f (x)
hm−1
2
xm−1
La riga i-esima del sistema lineare risulta dunque essere


..


.


..


.
i uj−1  
h

1
1
1
1
 uj  = f (xi ) hi−1 +hi 
+
−
0
.
.
.
0
0 . . . 0 − hi−1
hi−1
hi
hi
 


2
uj+1 
..


.
..
.
57
8.2. METODI SPETTRALI
e dunque uguale (per questo semplice problema modello) a quella della discretizzazione con differenze finite del secondo ordine. Pertanto, è naturale aspettarsi,
sotto opportune ipotesi di regolarità, che l’errore rispetto alla soluzione analitica
tenda a zero come h2 , h = maxj hj .
A questo punto si risolve il sistema lineare, dopo aver opportunamente modificato la matrice e il termine noto per imporre le condizioni al bordo di Dirichlet.
Nel caso di condizioni di Neumann omogenee, non è necessaria alcuna modifica.
Infatti, la forma debole del problema è
Z 1
¯1 Z 1
¯
′
′
′
f (x)φi (x)dx, i = 1, . . . , m
u (x)φi (x)dx =
−u (x)φi (x)¯ +
0
0
0
Il primo termine è naturalmente zero per 1 < i < m. Non considerarlo, cioè porlo
a zero, neanche per i = 1 (i = m) significa richiedere u′ (x1 ) = 0 (u′ (xm ) = 0) visto
che φ1 (x1 ) = 1 (φm (xm ) = 1). Nel caso di condizioni di Neumann non omogenee
(per esempio in 0), basta modificare la prima riga del sistema secondo l’equazione
Z 1
Z 1
f (x)φ1 (x)dx
u′ (x)φ′1 (x)dx = −u′ (0) +
0
0
8.1.2
Estensione al caso bidimensionale
Tutto quanto detto si estende, in particolare, al caso bidimensionale. Si deve usare
la formula di Green
Z
Z
Z
2
v(s)∇u(s) · ν(s)ds
∇ u(x)v(x)dx = − ∇u(x) · ∇v(x)dx +
Ω
∂Ω
Ω
ove ν(s) è il versore esterno normale a ∂Ω.
8.2
Metodi spettrali
Sia
u(x) =
X
uj φj (x)
j
L’indice algebrico di convergenza è il più grande k tale che
lim |uj |j k < +∞
j→∞
Se tale limite è finito per ogni k, allora la serie si dice convergere esponenzialmente
oppure spettralmente. Significa che |uj | decade più velocemente di ogni potenza
negativa di j. Parleremo di metodi spettrali quando useremo un’approssimazione
di una serie convergente spettralmente
u(x) ≈
m
X
j=1
ûj φj (x)
58
CAPITOLO 8. METODI VARIAZIONALI
per u(x).
Consideriamo un sistema {φj }j ortonormale rispetto al prodotto scalare
Z b
φj (x)φi (x)w(x)dx = δji
a
La formulazione di Galerkin di un problema ai limiti Lu = f , L operatore differenziale lineare, diventa
Z b
Z b
m
X
Lφj (x)φi (x)w(x)dx =
f (x)φi (x)w(x)dx
uj
a
j=1
a
Per l’approssimazione degli integrali si usano le corrispondenti formule di quadratura gaussiana a m punti, dando origine al sistema lineare
Ãm
!
m
m
X
X
X
uj
f (xn )φi (xn )wn
(8.6)
Lφj (xn )φi (xn )wn =
n=1
n=1
j=1
In tal caso si parla di metodi pseudospettrali. I coefficienti uj che si trovano risolvendo il sistema lineare si chiamano solitamente soluzione nello
P spazio spettrale.
Dati i coefficienti, si ricostruisce la soluzione nello spazio fisico j uj φj (x). La formulazione di Galerkin è particolarmente utile quando le funzioni {φj }j soddisfano
naturalmente le condizioni al bordo (condizioni al bordo periodiche, condizioni di
Dirichlet omogenee, vanishing boundary conditions) e/o il calcolo di Lφj (x) risulta
semplice.
Consideriamo ora un caso particolare di fondamentale importanza. Molte
proprietà risultato comuni anche agli altri metodi presudospettrali.
8.2.1
Trasformata di Fourier
Sia [a, b] un intervallo di R, m > 0 pari e fissato. Consideriamo, per ogni j ∈ Z,
φj (x) =
Allora,
Z
ei(j−1−m/2)2π(x−a)/(b−a)
√
b−a
b
φj (x)φk (x)dx = δjk
a
Infatti, se j = k allora φj (x)φk (x) = 1/(b − a), altrimenti
φj (x)φk (x) =
ei2π(j−k)(x−a)/(b−a)
b−a
e quindi
Z
a
b
φj (x)φk (x)dx =
Z
1
0
ei2π(j−k)y
(b − a)dy = 0 ,
b−a
(8.7)
59
8.2. METODI SPETTRALI
poiché l’integrale delle funzioni sin e cos in un intervallo multiplo del loro periodo è
nullo. La famiglia di funzioni {φj }j si dice ortonormale nell’intervallo [a, b] rispetto
al prodotto scalare
Z b
(φj , φk ) =
φj (x)φk (x)dx
a
Un risultato utile è il seguente
m
X
ei(n−1)2π(j−k)/m = mδjk ,
n=1
1 ≤ j, k ≤ m
(8.8)
È ovvio per j = k; altrimenti
m
X
e
i(n−1)2π(j−k)/m
=
m−1
X³
ei2π(j−k)/m
n=0
n=1
=
´n
=
1 − ei2π(j−k)
1 − cos(2π(j − k))
=
=0
i2π(j−k)/m
1−e
1 − ei2π(j−k)/m
poiché −m + 1 ≤ j − k ≤ m − 1.
8.2.2
Trasformata di Fourier discreta
Sia u una funzione da [a, b] a C periodica (u(a) = u(b)). Supponiamo che u si
possa scrivere (ciò è vero, per esempio, per funzioni di classe C 1 ) come
+∞
X
u(x) =
uj φj (x),
j=−∞
uj ∈ C
(8.9)
Fissato k ∈ Z, moltiplicando entrambi i membri per φk (x) e integrando nell’intervallo [a, b], usando (8.7) si ottiene
Z
b
u(x)φk (x)dx =
a
=
Z
b


a
∞
X
+∞
X
j=−∞
uj
j=−∞
Z

uj φj (x)φk (x) dx =
b
(8.10)
φj (x)φk (x)dx = uk
a
Dunque, abbiamo un’espressione esplicita per uj . Analogamente si vede che
Z
a
b
|u(x)|2 dx =
+∞
X
|uj |2
j=−∞
(identità di Parseval)
60
CAPITOLO 8. METODI VARIAZIONALI
La prima approssimazione da fare consiste nel troncare la serie infinita. Osserviamo
che, definito J = Z \ {1, 2, . . . , m},
¯
¯2
¯2
¯
¯
¯
Z b ¯X
Z b¯
m
X
¯
¯
¯
¯
¯
¯
¯u(x) −
uj φj (x)¯¯ dx =
uj φj (x)¯ dx =
¯
¯
a ¯ j∈J
a ¯
¯
¯
j=1

Ã
!
Z b X
X
X

|uk |2
uk φk (x) dx =
=
uj φj (x)
a
j∈J
k∈J
k∈J
Stimiamo adesso uk : si ha, per funzioni di classe C 1 ,
uk =
Z
a
b
µ
¶
b−a
u(b)
u(a)
√
u(x)φk (x)dx = −
+
−√
i(k − 1 − m/2)2π
b−a
b−a
Z b
b−a
u′ (x)φk (x)dx = O(k −1 )
+
i(k − 1 − m/2)2π a
Se anche u′ (a) = u′ (b) e u′ ∈ C 1 , allora, integrando ancora per parti, si ottiene
uk = O(k −2 ) e cosı̀ via. Se dunque u è infinitamente derivabile e periodica (cioè
tutte le derivate sono periodiche), allora uk decade più velocemente di ogni potenza
negativa di k.
La seconda approssimazione da fare è utilizzare una formula di quadratura per
il calcolo di uk . Riportiamo per comodità la formula di quadratura trapezoidale a
m + 1 nodi equispaziati xn = (b − a)yn + a, ove yn = (n − 1)/m, n = 1, . . . , m + 1
per funzioni periodiche:
!
Ã
Z b
m
m
X
b−a
b−aX
u(x)dx ≈
u(xn ) + f (xm+1 ) =
u(xn )
f (x1 ) + 2
2m
m
a
n=2
n=1
Usando la (8.8), abbiamo
mδjk =
m
X
ei(n−1)2π(j−k)/m =
m
X
ei(j−k)2πyn =
ei(j−k)2π(xn −a)/(b−a) =
n=1
n=1
n=1
m
X
m
X
ei(j−1−m/2)2π(xn −a)/(b−a) e−i(k−1−m/2)2π(xn −a)/(b−a)
√
√
= (b − a)
=
b
−
a
b
−
a
n=1
Z b
m
X
φj (x)φk (x)dx
φj (xn )φk (xn ) = m
= (b − a)
n=1
a
cioè la famiglia {φj }j , 1 ≤ j ≤ m, è ortonormale anche rispetto al prodotto scalare
discreto
m
b−aX
φj (xn )φk (xn )
(φj , φk )d =
m
n=1
61
8.2. METODI SPETTRALI
Questo significa che la formula di quadratura trapezoidale a m punti è esatta per
m−1
le funzioni {φj }j=−m+1
. Applicando la formula di quadratura ai coefficienti (8.10)
si ottiene
Z b
e−i(k−1−m/2)2π(x−a)/(b−a)
√
dx =
u(x)
uk =
b−a
a
Z 1
√
u((b − a)y + a)e−i(k−1)2πy eimπy dy ≈
= b−a
≈
√
b−a
m
0
m
X
n=1
¡
¢
u(xn )eimπyn e−i(k−1)2πyn = ûk
ove x = (b − a)y + a.
La funzione (serie troncata di Fourier)
ũ(x) =
m
X
m/2−1
ûj φj (x) =
j=1
X
m/2−1
=
X
ûk+1+m/2 φk+1+m/2 (x) =
k=−m/2
ûk+1+m/2
k=−m/2
eik2π(x−a)/(b−a)
√
b−a
è un polinomio trigonometrico che approssima f (x) ed è interpolante nei nodi xn .
Infatti, usando (8.8),
ũ(xn ) =
m
X
ûj φj (xn ) =
j=1
!
Ã√
m
¢ −i(j−1)2πy
ei(j−1−m/2)2π(xn −a)/(b−a)
b−a X¡
imπyk
k
√
u(xk )e
e
=
=
m
b−a
j=1
m
X
=
=
1
m
1
m
m
X
k=1
m
X
k=1
u(xk )eimπ(k−1)/m e−imπ(n−1)/m
m
X
e−i(j−1)2π(k−1)/m ei(j−1)2π(n−1)/m =
j=1
u(xk )ei(k−n)π
m
X
j=1
k=1
ei(j−1)2π(n−k)/m =
1
u(xn )m = u(xn ) .
m
Si può far vedere poi che
¯
¯2
¯
Z b¯
m
+∞
X
X
¯
¯
¯u(x) −
¯
|uk |2
ûj φj (x)¯ dx ≤ 2
¯
a ¯
¯
j=1
k∈J
La trasformazione
[u(x1 ), u(x2 ), . . . , u(xm )]T → [û1 , û2 , . . . , ûm ]T
62
CAPITOLO 8. METODI VARIAZIONALI
si chiama trasformata di Fourier discreta
√ di u e û1 , . . . , ûm coefficienti di Fourier
T
di u. Il vettore m · [û1 , û2 , . . . , ûm ] / b − a può essere scritto come prodotto
matrice-vettore F [u(x1 )eimπy1 , u(x2 )eimπy2 , . . . , u(xm )eimπym ]T , ove
F = (fjk ),
fjk = e−i(j−1)2πyk .
Alternativamente, si può usare la Fast Fourier Transform (FFT). Il comando fft
imπy1 , u(x )eimπy2 , . . . , u(x )eimπym ]T produce il vettore
applicato al vettore [u(x
m
2
√ 1 )e
T
m · [û1 , û2 , . . . , ûm ] / b − a, cosı̀ come il comando fftshift applicato al risultato
del comando fft applicato al vettore [u(x1 ), u(x2 ), . . . , u(xm )].
Dati dei coefficienti v̂j , j = 1, . . . , m, si può considerare la funzione (periodica)
m
X
v̂j φj (x)
j=1
La valutazione nei nodi xn , n = 1, . . . , m, porge
v̂ˆn =
m
X
v̂j φj (xn ) =
La trasformazione
v̂j
j=1
j=1
=√
m
X

m
ei(j−1−m/2)2π(xn −a)/(b−a)
√
=
b−a

1 X i(j−1)2πyn  −imπyn
m
.
e
v̂j e
b − a m j=1
[v̂1 , v̂2 , . . . , v̂m ]T → [v̂ˆ1 , v̂ˆ2 , . . . , v̂ˆm ]T
si chiama anti-trasformata di Fourier discreta. Se i v̂j sono i coefficienti di Fourier
di una funzione
di interpolazione comporta v̂ˆn = u(xn ).
√ v(x), la proprietà
imπy
ˆ
ˆ
1
Il vettore b − a · [v̂1 e
, v̂2 eimπy2 , . . . , v̂ˆm eimπym ]T /m può essere scritto co∗
me prodotto matrice-vettore F [v̂1 , v̂2 , . . . , v̂m ]T /m. Alternativamente,
il comando
√
imπy
ˆ2 eimπy2 , . . . , v̂ˆm eimπym ]/
ˆ
1 , v̂
ifft applicato al vettore [v̂1 , v̂2 , . . . , v̂m ] produce il vettore b − a·[v̂1 e
mentre, se applicato al risultato √
del comando ifftshift applicato al vettore
[v̂1 , v̂2 , . . . , v̂m ], produce il vettore b − a · [v̂ˆ1 , v̂ˆ2 , . . . , v̂ˆm ]/m.
Costi computazionali e stabilità
La Fast Fourier Transform di un vettore di lunghezza m ha costo O(m log m), mentre il prodotto matrice-vettore O(m2 ). Tali costi sono però asintotici e nascondono
i fattori costanti. Inoltre, GNU Octave può far uso di implementazioni ottimizzate
di algebra lineare (come, ad esempio, le librerie ATLAS). In pratica, dunque, esiste
un m0 sotto il quale conviene, dal punto di vista del costo computazionale, usare
il prodotto matrice-vettore e sopra il quale la FFT.
Per quanto riguarda l’accuratezza, in generale la FFT è più precisa del prodotto matrice vettore. Poiché la trasformata di Fourier discreta comporta l’uso
63
8.2. METODI SPETTRALI
di aritmetica complessa (anche se la funzione da trasformare è reale), la sequenza trasformata/anti-trasformata potrebbe introdurre una quantità immaginaria
spuria che può essere eliminata con il comando real.
Anche per la trasformata di Fourier vi possono essere problemi di stabilità
simili al fenomeno di Runge (qui chiamato fenomeno di Gibbs). Una tecnica per
“smussare” (in inglese “to smooth”) eventuali oscillazioni, consiste nel moltiplicare
opportunamente i coefficienti di Fourier ûj per opportuni termini σj che decadono
in j, per esempio
Ã
!
2|j − m/2 − 21 | − 1
σj = 1 −
, j = 1, . . . , m
(8.11)
m
Il risultato è che i coefficienti ûm/2 e ûm/2+1 sono pesati da σm/2 = σm/2+1 = 1,
i coefficienti fˆm/2−1 e fˆm/2+2 da σm/2−1 = σm/2+2 = m−2
m e cosı̀ via secondo
la formula (8.11). Questa scelta corrisponde alle medie di Cesáro. Infatti, si
sostituisce la serie troncata di Fourier
m
X
ũ(x) =
ûj φj (x)
j=1
Fk =
M
+k
X
fˆn φn (x)
k−1
X
fˆm+1+M φm+1+M (x),
k=M
m=−k
m=M +1−k
con la media delle troncate
m/2
X
k=1


m/2+k
X
j=m/2+1−k

ûj φj (x) /m
Si ricorda che se una serie è convergente, allora il limite delle medie delle sue
troncate è la somma della serie.
Valutazione di un polinomio trigonometrico
Supponiamo di conoscere i coefficienti ûj , j = 1, . . . , m e di voler valutare la
funzione
m
X
ûj φj (x) .
u(x) =
j=1
su un insieme di nodi target xk equispaziati, xk = (k − 1)/n, k = 1, . . . , n, n > m,
n pari. Si possono introdurre dei coefficienti fittizi Ûk
Ûk = 0
k = 1, . . . ,
Ûk = ûk− n−m
k=
2
Ûk = 0
n−m
2
n−m
n−m
+ 1, . . . , m −
2
2
n−m
k =m−
+ 1, . . . , n
2
64
CAPITOLO 8. METODI VARIAZIONALI
Si avrà
ˆk =
û
m
X
ûj φj (xk ) =
j=1
n
X
Ûj
j=1
ei(j−1−n/2)2π(xk −a)/(b−a)
√
=
b−a


n
1 X
n
Ûj ei(j−1)2πyk  e−inπyk
=√
b − a n j=1
Oppure si può costruire la matrice F relativa ai nodi (ciò funziona anche per nodi
non equispaziati). Infine, si può usare le NFFT.
8.3
Metodi di collocazione
Si assume comunque
u(x) =
m
X
uj φj (x)
j=1
ove {φj } è un sistema ortonormale rispetto ad un prodotto scalare, ma si impone
poi che l’equazione differenziale Lu = f sia soddisfatta in certi nodi xn . Si ha il
seguente risultato interessante:
Teorema 5. La soluzione del sistema lineare
m
X
uj Lφj (xn ) = f (xn ), n = 1, 2, . . . , m
(8.12)
j=1
ove gli {xn } sono i nodi della quadratura gaussiana relativa alla famiglia {φj } è
la stessa del problema di Galerkin
Z b
Z b
m
X
f (x)φi (x)w(x)dx
Lφj (x)φi (x)w(x)dx =
uj
j=1
a
a
quando si approssimino gli integrali con le formule gaussiane.
Dimostrazione. Per ogni i, 1 ≤ i ≤ m, da (8.12), si ha
m
X
uj Lφj (xn )φi (xn )wn = f (xn )φi (xn )wn ,
n = 1, 2, . . . , m
j=1
ove i {wn }n sono i pesi di quadratura gaussiana, da cui, sommando su n,


!
Ãm
m
m
m
X
X
X
X

Lφj (xn )φi (xn )wn =
uj Lφj (xn )φi (xn )wn  =
uj
n=1
j=1
j=1
=
m
X
n=1
f (xn )φi (xn )wn ,
i = 1, 2, . . . , m
n=1
che è precisamente la formulazione di Galerkin pseudospettrale (8.6).
65
8.3. METODI DI COLLOCAZIONE
8.3.1
Condizioni al bordo
Consideriamo il problema


 Lu(x) = f (x)
u(a) = α

 ′
u (b) = β
Con il metodo di collocazione, si ha
X
m


uj Lφj (xn ) = f (xn ),




j=1




m
X
uj φj (a) = α


j=1



m

X



uj φ′j (b) = β


n = 1, 2, . . . , m − 2
j=1
Anche in questo caso il metodo di collocazione può essere visto come un metodo
di Galerkin pseudospettrale: basta prendere come nodi di collocazione gli m − 2
nodi di quadratura gaussiana. Si ha poi

Ãm−2
! m−2
m
X
X
X




uj
f (xn )φi (xn )wn , i = 1, 2, . . . , m − 2
Lφj (xn )φi (xn )wn =




n=1
n=1
j=1



m
X
uj φj (a) = α


j=1




m

X



uj φ′j (b) = β


j=1
Alternativamente, si possono usare, come nodi di collocazione, quelli delle formule
di quadratura di Gauss–Lobatto (che contengono i nodi al bordo).
Collocazione Gauss–Lobatto–Chebyshev
I polinomi di Chebyshev sono definiti da
Tj (x) = cos(j arccos(x)),
e soddisfano
Z
1
−1



ππ
Tj (x)Ti (x)
√
dx =
2

1 − x2

0
−1 ≤ x ≤ 1
i=j=0
i = j 6= 0
i 6= j
66
CAPITOLO 8. METODI VARIAZIONALI
(lo si vede con il cambio di variabile x = cos θ e applicando le formule di Werner). I
nodi di Gauss–Lobatto sono xn = cos((n−1)π/(m−1)), n = 1, 2, . . . , m. Possiamo
allora definire la seguente famiglia di funzioni ortonormali
φ1 (x) =
r
1
T0 (x),
π
φj (x) =
r
2
Tj−1 (x),
π
j = 2, 3, . . . , m
Ricordando la formula di ricorrenza tra polinomi di Chebyshev, possiamo scrivere
r
√
1
2
φ1 (x) =
, φ2 (x) =
x, φ3 (x) = 2xφ2 (x) − 2φ1 (x),
π
π
φj+1 (x) = 2xφj (x) − φj−1 (x), j = 3, 4, . . . , m − 1
r
Da qui, possiamo calcolare anche la derivata prima e seconda delle funzioni:
r
2
, φ′3 (x) = 2φ2 (x) + 2xφ′2 (x),
π
φ′j+1 (x) = 2φj (x) + 2xφ′j (x) − φ′j−1 (x), j = 3, 4, . . . , m − 1
φ′1 (x)
= 0,
φ′′1 (x) = 0,
φ′2 (x)
=
φ′′2 (x) = 0,
φ′′3 (x) = 4φ′2 (x),
φ′′j+1 (x) = 4φ′j (x) + 2xφ′′j (x) − φ′′j−1 (x),
Conviene calcolare le matrici

φ1 (x1 )
 φ2 (x1 )

T= .
 ..
φ1 (x2 )
φ2 (x2 )
..
.
φm (x1 ) φm (x2 )
 ′
φ1 (x1 ) φ′1 (x2 )
 φ′ (x1 ) φ′ (x2 )
2
 2
T′ =  .
..
 ..
.
φ′m (x1 ) φ′m (x2 )
 ′′
φ1 (x1 ) φ′′1 (x2 )
 φ′′ (x1 ) φ′′ (x2 )
2
 2
T′′ =  .
..
 ..
.
...
...
j = 3, 4, . . . , m − 1

φ1 (xm )
φ2 (xm ) 


..

.
...
. . . φm (xm )

. . . φ′1 (xm )
. . . φ′2 (xm ) 


..

...
.
. . . φ′m (xm )

. . . φ′′1 (xm )
. . . φ′′2 (xm ) 


..

.
...
φ′′m (x1 ) φ′′m (x2 ) . . . φ′′m (xm )
Consideriamo, a titolo di esempio, il seguente problema modello

′′

 −u (x) + q(x)u(x) = f (x)
u(−1) = α


u′ (1) = β
67
8.3. METODI DI COLLOCAZIONE
Il sistema lineare da risolvere per il metodo di collocazione Gauss–Chebyshev–
Lobatto (per il momento senza tener conto delle condizioni al bordo) è

   


f (x1 )
u1
q(x1 )
   u2   f (x2 ) 


q(x2 )

 T   
 ′′ T 
−T
+
 T   ..  =  .. 


..







.
. 
.
f (xm )
um
q(xm )
Per imporre le condizioni al bordo, si sostituisce la prima riga della matrice con
la prima riga di TT e il primo elemento del termine noto con α. Poi, l’ultima riga
della matrice con l’ultima riga di T′ T e l’ultimo elemento del termine noto con
β. Una volta noti i coefficienti uj , si ricostruisce la soluzione nello spazio fisico
tramite
 


u1
u(x1 )
 u2 
 u(x2 ) 
 


 ..  = TT  .. 
 . 
 . 
u(xm )
um
Capitolo 9
Esercizi
1. Si risolva il problema ai limiti
 ′′

 u (x) = u(x) + x,
u(0) = 0


u(1) = 0
x ∈ (0, 1)
(9.1)
usando il metodo delle differenze finite del secondo ordine. Sapendo che la
soluzione esatta è u(x) = (ex − e−x )/(e − e−1 ) − x, si mostri inoltre l’ordine
del metodo mediante un grafico logaritmico-logaritmico dell’errore in norma
infinito. Si risolva infine lo stesso problema con il metodo di shooting.
2. Si risolva il problema ai limiti
 ′′
′

 u (x) + u (x) + u(x) − cos(x) = 0,
u(0) = 0


u(π/2) = 1
x ∈ (0, π/2)
usando il metodo delle differenze finite del secondo ordine. Si mostri inoltre
l’ordine del metodo mediante un grafico logaritmico-logaritmico dell’errore
in norma infinito rispetto ad una soluzione di riferimento.
3. Si risolva il problema ai limiti
 ′′
′

 u (x) + u (x) + u(x) − cos(x) = 0,
u′ (0) = 1


u(π/2) = 1
x ∈ (0, π/2)
usando il metodo delle differenze finite del secondo ordine. Si mostri inoltre
l’ordine del metodo mediante un grafico logaritmico-logaritmico dell’errore
in norma infinito rispetto ad una soluzione di riferimento.
68
69
4. Si risolva il problema ai limiti
 ′′

 u (x) = cos(u(x)),
u(0) = 0


u(1) = 1
x ∈ (0, 1)
usando il metodo delle differenze finite del secondo ordine.
5. Si risolva il problema ai limiti
µ
¶

d
d



 − dx (1 + x) dx u(x) = 1,
u(0) = 0




u(1) = 0
x ∈ (0, 1)
Si mostri inoltre l’ordine del metodo mediante un grafico logaritmico-logaritmico
dell’errore in norma infinito rispetto alla soluzione esatta.
6. Si risolva il problema ai limiti
 ′′
′

 u (x) = 20u (x) + u(x),
u(0) = 0


u(1) = 1
x ∈ (0, 1)
Visto l’andamento della soluzione, si implementi uno schema di differenza
finite su nodi non equispaziati secondo una distribuzione di tipo coseno.
7. Si ricavi la relazione di ricorrenza dei polinomi ortonormali nell’intervallo
2 2
[−∞, ∞] rispetto alla funzione peso w(x) = e−α x
8. Noti gli zeri dei polinomi di Legendre e i pesi di quadratura della rispettiva formula gaussiana, si ricavino i nodi e i pesi di una formula gaussiana
nell’intervallo [a, b] rispetto al peso w(x) = 1.
9. Si risolva il problema ai limiti (9.1) usando il metodo di collocazione con
polinomi di Legendre. Gli N nodi di collocazione in [a, b] e la valutazione dei polinomi di Legendre e delle loro derivate sono dati dalla function
[L,x,L1,L2] = legendrepolynomials(N,a,b). Si mostri inoltre l’ordine
del metodo mediante un grafico logaritmico-logaritmico dell’errore in norma
infinito.
Parte III
PDEs
(Equazioni alle derivate
parziali)
70
Capitolo 10
Equazione del calore
10.1
Equazione del calore con dati iniziali e
condizioni ai limiti
Consideriamo la seguente equazione alle derivate parziali

∂2u
∂u



(t,
x)
=
(t, x), t > 0, x ∈ (0, L)
 ∂t
∂x2
u(a, t) = u(b, t) = 0,
t > 0 (condizioni ai bordi)



 u(0, x) = u (x),
x ∈ (0, L) (condizioni iniziali)
0
(10.1)
Supponiamo che u0 (x) verifichi le condizioni di compatibilità u0 (a) = u0 (b) = 0.
Tale equazione rappresenta, per esempio, l’andamento della temperatura u su una
barra di lunghezza L, i cui estremi sono tenuti a temperatura zero, e con una
iniziale distribuzione di temperatura u0 (x).
10.1.1
Esistenza di una soluzione
Cerchiamo una soluzione a variabili separabili
u(t, x) = ψ(t)φ(x)
Inserendo tale rappresentazione in (10.1), si deduce
ψ ′ (t)φ(x) = ψ(t)φ′′ (x),
t > 0, x ∈ (0, L)
da cui
ψ ′ (t)
= −K (costante) ⇒ ψ(t) = Ae−Kt
ψ(t)
Per quanto riguarda φ(x), l’unica possibilità non banale, tenendo conto anche delle
condizioni ai bordi, è che K = λ2 > 0, λ > 0. Allora
φ′′ (x)
= −λ2 ⇒ φ(x) = B sin(λx) + C cos(λx)
φ(x)
71
72
CAPITOLO 10. EQUAZIONE DEL CALORE
Le condizioni ai bordi impongono C = 0 e λ = kπ/L, k numero naturale non nullo.
Pertanto, la funzione
µ 2 2 ¶
µ
¶
k π
kπ
uk (t, x) = exp − 2 t sin
x
L
L
è soluzione dell’equazione del calore (e soddisfa le condizioni ai bordi) per ogni k.
Quindi, la seguente combinazione lineare infinita
u(t, x) =
∞
X
Ak uk (t, x)
k=1
è soluzione formale dell’equazione del calore. Per quanto riguarda la condizione
iniziale, si deve imporre
µ
¶
∞
X
kπ
Ak sin
u0 (x) = u(0, x) =
x
L
k=1
Possiamo prolungare per antisimmetria la funzione u0 (x) nell’intervallo [−L, L].
Sotto opportune ipotesi, la sua serie di Fourier
ū0 (x) =
+∞
X
u0 j φj (x)
j=−∞
converge. Si ha dunque, per parità,
√ Z
Z L
−i
−i 2 L
um/2+1+k = √
ū0 (x) sin(2πk(x + L)/(2L))dx = √
ū0 (x) sin(kπx/L + kπ)dx
2L −L
L 0
e
um/2+1−k
i
=√
2L
Z
√ Z
i 2 L
ū0 (x) sin(2πk(x + L)/(2L))dx = √
ū0 (x) sin(kπx/L + kπ)dx
L 0
−L
L
da cui
ū0 (x) =
=
+∞
X
k=−∞
∞
X
um/2+1+k φm/2+1+k (x) =
um/2+1+k
k=−∞
∞ √
X
cos(kπx/L + kπ) + i sin(kπx/L + kπ)
√
=
2L
2
√ um/2+1+k i sin(kπx/L + kπ) =
L
k=1
µ
¶ ¸
µ
¶
·
∞
X 2Z L
kπ
kπ
u0 (x) sin
x dx sin
x
=
L 0
L
L
=
k=1
10.1. EQUAZIONE DEL CALORE CON DATI INIZIALI E CONDIZIONI AI LIMITI73
Si potrebbe mostrare adesso che
u(t, x) =
∞ · Z
X
2
L
k=1
L
u0 (x) sin
0
¶ ¸
µ 2 2 ¶
µ
¶
kπ
kπ
k π
x dx exp − 2 t sin
x
L
L
L
µ
è soluzione di (10.1). Dalla presenza del termine esponenziale negativo nel tempo,
si deduce che la soluzione tende a zero per t → +∞. L’equazione del calore rappresenta il modello dei fenomeni di diffusione. La diffusione è il processo mediante il
quale la materia (o l’energia) è trasportata da una parte di un sistema ad un’altra
come risultato di moti molecolari random. Notiamo come la funzione
U (t, x) = u(−t, x)
per cui limt→+∞ U (t, x) = ∞, soddisfi l’equazione differenziale
∂U
∂2U
(t, x) = − 2 (t, x)
∂t
∂x
con le stesse condizioni ai bordi e iniziali dell’equazione del calore.
10.1.2
Unicità della soluzione
Introduciamo la seguente quantità (energia)
E(t) =
Z
L
0
1 2
u (t, x)dx
2
Si ha
dE
=
dt
Z
L
0
¸
·
Z L
Z L 2
∂ 1 2
∂u
∂ u
u dx =
u 2 dx
u (t, x) dx =
∂t 2
∂t
∂x
0
0
Integrando per parti e tenendo conto delle condizioni ai bordi, si ha
dE
=−
dt
Z
0
Lµ
∂u
∂x
¶
dx ≤ 0
Per dimostrare l’unicità, consideriamo come al solito il problema omogeneo (corrispondente a (10.1) con u0 ≡ 0. Per tale problema E0 (0) = 0 e quindi 0 ≤ E0 (t) ≤
E0 (0) da cui E0 (t) = 0 per ogni t. Quindi u(t, x) ≡ 0 è l’unica soluzione del problema omogeneo. Dunque, se u1 (t, x) e u2 (t, x) fossero due soluzioni del problema
(10.1), allora u1 (t, x) − u2 (t, x) sarebbe soluzione del problema omogeneo e quindi
u1 (t, x) ≡ u2 (t, x).
74
10.2
CAPITOLO 10. EQUAZIONE DEL CALORE
Metodo delle linee
Il metodo delle linee per la risoluzione di problemi del tipo

∂u
∂2u



(t,
x)
=
(t, x) + f (u(t, x)) + g(t, x), t > 0, x ∈ (a, b)
 ∂t
∂x2
 + condizioni ai bordi


 + condizioni iniziali
(10.2)
ove il termine f (u(t, x)) si chiama reazione e il termine g(t, x) sorgente, prevede
di discretizzare gli operatori differenziali spaziali con uno dei metodi visti per i
problemi con valori ai bordi e poi risolvere il sistema di ODEs che ne risulta con
un metodo per problemi ai valori iniziali visti. Vediamo qualche esempio.
10.2.1
Differenze finite
Trascurando per il momento le condizioni ai bordi e usando differenze finite centrate del secondo ordine a passo costante h



 ′
 
 

−2
1
0
...
0 
y1 (t)
y1 (t)
f (y1 (t))
g1 (t)


.
..   y (t)   f (y (t))   g (t) 
 y ′ (t) 
 1 −2
1
...
2

 2
 
 2
  2


1





 



..
.
.
..
.
.
.
.
.
.
.
.
=
+
+










.
.
.
.
.
0 
.
 h2  0
 ′ .
 
 

 ym−1 (t) f (ym−1 (t)) gm−1 (t)
 ..
ym−1 (t)
..


.
1
−2
1
.
′ (t)
ym
ym (t)
f (ym (t))
gm (t)
0
... 0
1 −2
ove yj (t) ≈ y(t, xj ) o, in maniera compatta,
y ′ (t) = Ay(t) + f (y(t)) + g(t)
(10.3)
(con l’ovvia definizione dei simboli). A questo punto, si sceglie il metodo di integrazione temporale (θ-metodo, Runge–Kutta, multistep). Si tenga presente che
il problema (10.3), che si dice semidiscretizzato, è solitamente un problema stiff.
Infine, si sistemano le condizioni ai bordi.
Condizioni al bordo di Dirichlet
Vediamo come imporre una condizione di Dirichlet costante in x1 = a. La prima
equazione del sistema di ODEs riguarda y1 (t) ≈ y(t, x1 ) e sarà della forma y1′ (t) =
. . . Se si vuole imporre una condizione di Dirichlet costante, è sufficiente modifcare
opportunamente la prima del termine di destra del sistema di ODEs in modo da
ottenere l’equazione y1′ (t) = 0.
Nel caso si volesse imporre una condizione di Dirichlet dipendente dal tempo,
per esempio y(t, a) = α(t), si deve ancora modificare la prima riga del sistema differenziale, ma tale modifica dipende dal metodo scelto per l’integrazione temporale.
75
10.2. METODO DELLE LINEE
Per esempio, avendo scelto il metodo di Eulero implicito, si ha
(I − kA)y n+1 − kf (y n+1 ) = y n + kg(tn+1 )
ove k è il passo di integrazione temporale, cioè si deve risolvere il sistema non
lineare in y n+1
B(y n+1 )y n+1 = b
A questo punto, si può sostituire la prima riga di B(y n+1 ) con il primo vettore
della base canonica e la prima riga di b a α(tn+1 ).
Per il metodo di Eulero esponenziale, invece, si ha
y n+1 = exp(kA)y n + kϕ1 (kA)(f (y n ) + g(tn )) = exp(kA)y n + kϕ1 (kA)b
Se la prima riga di A viene messa a zero, la prima riga di exp(kA) e ϕ1 (kA) è
il primo vettore della base canonica e dunque basta porre il primo elemento di b
uguale a (α(tn+1 ) − α(tn ))/k.
Condizioni al bordo di Neumann
Per quanto riguarda una condizione di Neumann omogenea, per esempio in x = b,
si può pensare di introdurre la variabile fittizia ym+1 (t) ≈ u(t, xm+1 ), xm+1 = b+h
2
e imporre che ym+1 (t) = ym−1 (t). L’approssimazione da usare per ∂∂xu2 (t, b) diventa
dunque
∂2u
u(t, xm+1 ) − 2u(t, xm ) + u(t, xm−1 )
(t, b) ≈
=
∂x2
h2
ym+1 (t) − 2ym (t) + ym−1 (t)
2ym−1 (t) − 2ym (t)
=
=
2
h
h2
In maniera analoga si possono trattare condizioni di Neumann non omogenee.
10.2.2
Elementi finiti
Nel caso di discretizzazione spaziale con elementi finiti lineari, la discretizzazione
del problema (10.2) porta al sistema di ODEs
P y ′ (t) = Ay(t) + f (y(t)) + g(t)
(10.4)
ove A è (l’opposta de) la stiffness matrix e P la mass matrix, definita da,
Z xj+1
hj−1 + hj
φj (x)φj (x)dx =
pjj =
3
xj−1
Z xj+1
(10.5a)
hj
pj j+1 = pj+1 j =
φj (x)φj+1 (x)dx =
6
xj
76
CAPITOLO 10. EQUAZIONE DEL CALORE
mentre, per j = 1 e j = m,
Z x2
h1
p11 =
φ1 (x)φ1 (x)dx =
3
Zx1x2
h1
φ1 (x)φ2 (x)dx =
p12 =
6
x1
Z xm
hm−1
pm−1 m = pm m−1 =
φm (x)φm−1 (x)dx =
6
xm−1
Z xm
hm−1
pmm =
φm (x)φm (x)dx =
3
xm−1
Poi, per 1 < i < m,


Z xi+1
m
X
fi =
f
uj φj (x) φi (x)dx =
xi−1
=
Z
xi
xi−1
(10.5b)
j=1


Z
m
X
f
uj φj (x) φi (x)dx +
j=1
xi+1
xi
hi
hi−1 + hi
hi−1
+ f (ui ) = f (ui )
≈ f (ui )
2
2
2
Z xi+1
hi−1 + hi
gi =
g(t, x)φi (x)dx ≈ g(t, xi )
2
xi−1


m
X
f
uj φj (x) φi (x)dx ≈
j=1
mentre per i = 1 e i = m
hm−1
h1
,
fm = f (um )
2
2
hm−1
h1
g1 = g(t, x1 ) , gm = g(t, xm )
2
2
Le condizioni di Dirichlet omogenee per un nodo xi si impongono sostituendo la
riga corrispondente di P con zeri e 1 in diagonale e la riga corrispondente del
termine di destra del sistema di ODEs con zeri.
Usando un metodo esplicito per la risoluzione del sistema differenziale (10.4),
è necessaria l’inversione della matrice di massa. Per tale motivo, si può ricorrere
alla tecnica del mass lumping che consiste nel rendere diagonale la matrice P
sostituendo ogni sua riga con una riga di zeri e la somma degli elementi originali in
diagonale. Tale modifica è equivalente all’approssimazione degli integrali in (10.5)
mediante la formula dei trapezi e dunque non riduce l’accuratezza del metodo.
(−1)
Infatti, la matrice PL A (PL la matrice di massa con lumping) risulta uguale
alla matrice che si ottiene discretizzando con differenze finite centrate del secondo
ordine.
Usando invece un metodo implicito per la risoluzione del sistema differenziale (10.4), non è necessaria la tecnica del mass lumping: semplicemente, si devono
risolvere sistemi lineari in cui la matrice identità è sostituita dalla matrice di massa.
f1 = f (u1 )
10.3. ESERCIZI
10.3
Esercizi
1. Si risolva la PDE

∂u
∂2u

2 t


 ∂t (t, x) = ∂x2 (t, x) + (1 + π )e sin(πx), t > 0, x ∈ (0, 1)
u(t, 0) = u(t, 1) = 0
t>0



 u(0, x) = sin(πx)
x ∈ (0, 1)
usando il metodo delle linee. Si mostrino gli ordini spaziali e temporali.
77
Capitolo 11
Equazione del trasporto
11.1
Linee caratteristiche
Consideriamo la seguente equazione del trasporto

 ∂u (t, x) + c ∂u (t, x) = 0, t > 0, x ∈ R
∂t
∂x

u(0, x) = u0 (x),
x∈R
(11.1)
È facile verificare che la soluzione esatta è u(t, x) = u0 (x − ct) (da cui il nome
dell’equazione). Dunque, si può definire la linea caratteristica x(t) = x0 + ct la
soluzione è costante e vale u0 (x0 ). Ovviamente esiste una linea caratteristica per
ogni x0 ∈ R. Per un problema più generale

 ∂u (t, x) + c(x) ∂u (t, x) = f (x), t > 0, x ∈ R
∂t
∂x

u(0, x) = u0 (x),
x∈R
possiamo definire le linee caratteristiche come quelle curve parametriche x(t) per
cui la soluzione u(t, x(t)) soddisfa un’equazione differenziale ordinaria, cioè

∂u
∂u
d


u(t,
x(t))
=
(t,
x(t))
+
(t, x(t))x′ (t) = f (x(t)), t > 0


dt
∂t
∂x


u(0, x0 ) = u0 (x0 )


x′ (t) = c(x(t)),
t>0




x(0) = x0
Si ha dunque
u(t, x(t)) = u0 (x0 ) +
Z
t
f (x(s))ds
0
e dunque la soluzione u(t̄, x̄) dipende solo dai valori di u0 in x(0) e di f in x(t),
0 ≤ t ≤ t̄, ove
(
x′ (t) = c(x(t)), 0 ≤ t < t̄
x(t̄) = x̄
78
79
11.1. LINEE CARATTERISTICHE
L’insieme {(t, x(t)), 0 ≤ t ≤ t̄} è chiamato dominio di dipendenza continuo di
u(t̄, x̄). Tornando al caso dell’equazione del trasporto, il dominio di dipendenza
continuo di u(t̄, x̄) è la retta t = 1c (x(t) − x(0)) nel piano tOx.
c>0
t
t = 1c (x − x̄) + t̄
t̄
a
x̄
b
x
Figura 11.1: Linee caratteristiche per l’equazione del trasporto, c > 0
È ovviamente più fisico considerare un dominio limitato a < x < b. Nel caso
(11.1) con c > 0, si vede dalla Figura 11.1 allora che ha senso (ed è necessario)
prescrivere un’unica condizione al bordo in x = a. Tale punto si chiama punto di
inflow mentre il punto x = b è detto di outflow. L’equazione di trasporto su un
dominio limitato si scrive allora

∂u
∂u


 ∂t (t, x) + c ∂x (t, x) = 0, t > 0, x ∈ (a, b), c > 0
u(0, x) = u0 (x),



u(t, a) = f (t),
x ∈ [a, b]
t>0
con u0 (a) = f (0). La soluzione analitica è u(t, x) = ũ0 (x − ct), ove
(
u0 (x)
a≤x≤b
¡ a−x ¢
ũ0 (x) =
x<a
f c
(11.2)
80
11.2
CAPITOLO 11. EQUAZIONE DEL TRASPORTO
Differenze finite
Introduciamo la seguente discretizzazione per l’equazione di trasporto (11.2)
xj = a + jh, j = 1, 2, . . . , m, h = (b − a)/m
tn = nk, n = 0, 1, . . .
un j ≈ u(tn , xj )
un 0 = f (tn )
11.2.1
Eulero esplicito/upwind
Se discretizziamo con differenze finite in avanti nel tempo e all’indietro nello spazio,
otteniamo
un+1 j − un j
un j − un j−1
+c
=0
k
h
da cui
un+1 j = (1 − cλ)un j + cλun j−1 ,
λ=
k
h
(11.3)
Introducendo la notazione un = [un 1 , un 2 , . . . , un m ]T , si ha
un+1 = U un + cλfn
(11.4)
ove

1 − cλ
0
...
 cλ
1
−
cλ
0


.
.
..
..
U =
 0
 ..
..
..
 .
.
.
0
...
0

...
0
...
0 

.. 
..
.
. 
,

..
.
0 
cλ 1 − cλ


un 0
 0 


 .. 

fn =  . 

 .. 
 . 
0
Consistenza
Si ha
1
[u(tn+1 , xj ) − (1 − cλ)u(tn , xj ) − cλu(tn , xj−1 )] =
k"
∂u
ck
1
= u(tn , xj ) + k (tn , xj ) + O(k 2 ) − u(tn , xj ) + u(tn , xj )+
k
∂t
h
µ
¶#
ck
∂u
−
u(tn , xj ) −
(tn , xj )h + O(h2 )
= O(k + h) = τn+1 j
h
∂x
81
11.2. DIFFERENZE FINITE
Stabilità
Consideriamo due schemi (11.4) perturbati
ũn+1 = U ũn + cλfn
ûn+1 = U ûn + cλfn
Allora
(ũn+1 − ûn+1 ) = U (ũn − ûn ) = U n+1 (ũ0 − û0 )
da cui
kũn+1 − ûn+1 k ≤ kU kn+1 kũ0 − û0 k
Dunque, se kU k1 = kU k∞ = |1 − cλ| + |cλ| ≤ 1, cioè
c
k
≤1
h
(11.5)
si ha stabilità del metodo.
Convergenza
Definiamo
τn+1 = [τn+1 1, , τn+1 2 , . . . , τn+1 m ]T
u(tn+1 ) = [u(tn , x1 ), u(tn , x2 ), . . . , u(tn , xm )]T
Allora, dal risultato di consistenza,
u(tn+1 ) = U u(tn ) + λfn + kτn+1
e dunque, definendo
en+1 = un+1 − u(tn+1 )
si ha
en+1 = U en − kτn+1
da cui, se vale (11.5),
ken+1 k∞ ≤ ken k∞ + kkτn+1 k∞ ≤ ke0 k∞ + tn+1 O(k + h)
11.2.2
Condizione CFL
Analizziamo, alla luce delle linee caratteristiche, il significato della condizione (11.5).
Essa significa chiedere che il dominio di dipendenza discreto di un+1 j (rappresentato dai pallini neri in Figura 11.2) “contenga” il dominio di dipenda continuo di
u(tn+1 , xj ). La condizione, nel caso generale,
|c|h
≤1
k
82
CAPITOLO 11. EQUAZIONE DEL TRASPORTO
c>0
t
t = 1c (x − xj ) + tn+1
t = hk (x − xj ) + tn+1
tn+1
tn
tn−1
a
xj−2
xj−1
b
xj
x
Figura 11.2: Dominio di dipendenza discreto dello schema (11.3), c > 0.
c>0
t
t = 1c (x − xj ) + tn+1
tn+1
tn
tn−1
a
xj
xj+1
xj+2
b
x
Figura 11.3: Dominio di dipendenza discreto dello schema (11.6), c > 0.
chiama CFL condition o condizione di Courant–Friedrichs–Lewy ed è una condizione necessaria (ma non sufficiente) per la stabilità. Per esempio, discretizzando con
differenze finite in avanti sia nel tempo che nello spazio l’equazione del trasporto
83
11.2. DIFFERENZE FINITE
con c > 0, otteniamo
un+1 j = (1 + cλ)un j − cλun j+1
(11.6)
da cui il dominio di dipendenza discreto in Figura 11.3. Chiaramente, in tal caso
la soluzione numerica non può convergere alla soluzione analitica. Anche quando il
dominio di dipendenza discreto contiene quello continuo si possono avere problemi
di stabilità: se consideriamo infatti lo schema che si ottiene usando differenze finite
in avanti nel tempo e centrate nello spazio, si ottiene
un+1 = Cun + cλfn
ove

1
 λ
c
 2

C=0

 ..
 .
0
−c λ2
1
..
.
..
.
...
0
...
−c λ2
..
.
...
..
.
c λ2
0
(11.7)
0
..
.





0 


λ 
1
−c 2
cλ 1 − cλ
(l’ultima riga è stata corretta con differenze finite all’indietro). Si vede dunque
che kCk1 = kCk∞ = 1 + cλ > 1.
11.2.3
Equazioni modificate
Consideriamo l’approssimazione mediante differenze finite all’indietro del termine
∂u
∂x (t, xj ): si ha
un j − un j−1
∂u
(tn , xj ) + O(h) =
=
∂x
h
un j+1 − un j−1 h un j+1 − 2un j + un j−1
−
=
=
2h
2
h2
∂u
h ∂2u
=
(tn , xj ) −
(tn , xj ) + O(h2 )
∂x
2 ∂x2
(11.8)
Dunque essa corrisponde ad una approssimazione del secondo ordine in h di un’equazione di trasporto-diffusione, con coefficiente di diffusione proporzionale ad h.
Ciò giustifica l’andamento dissipativo della soluzione dell’equazione del trasporto
mediante lo schema upwind. In generale, la presenza di termini dissipativi stabilizza la discretizzazione, perché anche gli errori vengono dissipati. Analogamente,
discretizzando mediante differenze finite centrate il termine ∂u
∂x (t, xj ), si ottiene
una approssimazione al quarto ordine in h di
∂u
h2 ∂ 3 u
u(t, x) +
(t, x)
∂x
6 ∂x3
Ciò giustifica l’andamento dispersivo della soluzione del trasporto mediante lo
schema Eulero/centrato (11.7).
84
CAPITOLO 11. EQUAZIONE DEL TRASPORTO
11.2.4
Formulazione di flusso
Un qualsiasi metodo alle differenze finite può essere scritto nella forma
un+1 j = un j − λ(Fj+1/2 − Fj−1/2 )
(11.9)
con Fj+1/2 = F (un j , un j+1 ), ove F (·, ·) è detto flusso numerico ed è, in generale,
un’approssimazione del flusso esatto cu(tn , xj+1/2 ) associato all’equazione (11.2).
La (11.9) può essere interpretata come una discretizzazione in xj e tn dell’equazione
che deriva da
¤
1£
d
ūj (t) +
cu(t, xj+1/2 ) − cu(t, xj−1/2 ) = 0
dt
h
∂u
∂u
(t, x) + c (t, x) = 0
∂t
∂x
integrando nello spazio tra xj−1/2 e xj+1/2 e ponendo
Z
1 xj+1/2
u(t, x)dx
ūj (t) =
h xj−1/2
La scelta del flusso numerico caratterizza il metodo. Citiamo i seguenti metodi
espliciti, indicando per ognuno l’ordine dell’errore di troncamento:
Eulero in avanti/centrato (O(k + h2 )):
1
Fj+1/2 = c(un j+1 + un j )
2
Equivale al metodo (11.7).
Lax–Friedrichs (O(h2 /k + k + h2 )):
1
Fj+1/2 = [c(un j+1 + un j ) − λ−1 (un j+1 − un j )]
2
È una modifica di (11.7), ove si sostituisce un j con la media (un j+1 +
un j−1 )/2. Ciò rende il metodo stabile, in quanto la norma della matrice
associata vale 1 (qualora la CFL condition sia soddisfatta).
Lax–Wendroff (O(k 2 + h2 + h2 k)):
1
Fj+1/2 = [c(un j+1 + un j ) − λc2 (un j+1 − un j )]
2
upwind (O(k + h)):
1
Fj+1/2 = [c(un j+1 + un j ) − |c|(un j+1 − un j )]
2
Equivale a (11.3).
11.2. DIFFERENZE FINITE
11.2.5
85
Verifica dell’ordine
Per verificare l’ordine di un metodo di ordine p nel tempo e q nello spazio, in cui
non sia possibile disaccoppiare il passo temporale k e il passo spaziale h, è possibile
scegliere un passo di discretizzazione temporale k proporzionale alla potenza q/p
del passo di discretizzazione spaziale. In questo modo, l’errore globale risulta essere
q
p
O(k p + hq ) = O(h p + hq ) = O(hq ).
11.2.6
Estensione
I metodi descritti sopra si estendono, senza modifiche, al caso c < 0, ove però
le condizioni al bordo in (11.1) vanno imposte in x = b. Si estendono inoltre a
problemi a coefficienti non costanti in forma conservativa
∂u
∂(c(t, x)u(t, x))
(t, x) +
=0,
∂t
∂x
ove Fj+1/2 = F (xj+1/2 , un j , un j+1 ). Se c(t, x) > 0, le condizioni al bordo di inflow
sono in x = a, se c(t, x) < 0 le condizioni al bordo di inflow sono in x = b. Negli
altri casi, essendo c(t, x) la velocità di propagazione, si possono originare shocks.
11.2.7
Condizioni al bordo
Nel punto di outflow non è possibile usare una discretizzazione della derivata
che utilizzi punti esterni al dominio. Si possono allora usare discretizzazioni non
simmetriche, eventualmente con ordine di approssimazione minore (vedi § 6.3).
Oppure, si possono estrapolare i valori mancanti.
11.2.8
Problemi di convezione-diffusione
I problemi di convezione-diffusione

∂u
∂2u
∂u



(t,
x)
+
c
(t,
x)
=
d
(t, x), t > 0, x ∈ (a, b), d > 0
 ∂t
∂x
∂x2
u(0, x) = u0 (x)



 + condizioni ai bordi
si risolvono, generalmente, con il metodo delle linee. Per quanto riguarda la discretizzazione spaziale, allora, si può usare il metodo delle differenze finite, degli
elementi finiti o metodi spettrali. Per quanto riguarda i primi due, occorre però
un passo di discretizzazione spaziale h sufficientemente piccolo. Nel caso di differenze finite centrate del secondo ordine, omettendo la dipendenza temporale, si ha
86
CAPITOLO 11. EQUAZIONE DEL TRASPORTO
infatti (per c > 0)
c
uj+1 − uj−1
uj+1 − 2uj + uj−1
∂2u
∂u
(t, xj ) − d 2 (t, xj ) ≈ c
−d
=
2
∂x
∂x
2h
µ
¶ h
uj − uj−1
uj+1 − 2uj + uj−1
ch
=c
+
−d
h
2
h2
Se il coefficiente ch/2 − d fosse non negativo, lo schema corrisponderebbe alla
discretizzazione di un’equazione senza diffusione o addirittura con un termine
antidiffusivo. Pertanto, si richiede, nel caso generale,
|c|h
<1
2d
La quantità |c|h/(2d) si chiama numero di Peclét di griglia.
Le condizioni ai bordi, di Dirichlet o di Neumann, si impongono come nel caso
dell’equazione del calore.
11.3
Esercizi
1. Data l’approssimazione del secondo ordine della derivata y ′ (xN −1 ) ≈ (yN −
yN −2 )/(2h), si ricavi per estrapolazione lineare il valore yN da yN −2 e yN −1
e lo si inserisca nell’approssimazione. Determinare l’ordine dell’approssimazione risultante.
2. Data l’approssimazione del secondo ordine della derivata seconda y ′′ (xN −1 ) ≈
(yN − 2yN −1 + yN −2 )/(h2 ), si ricavi per estrapolazione quadratica il valore
yN da yN −3 , yN −2 e yN −1 e lo si inserisca nell’approssimazione. Determinare
l’ordine dell’approssimazione risultante.
3. Data l’approssimazione del secondo ordine della derivata y ′ (xN −1 ) ≈ (yN −
yN −2 )/(2h), si ricavi per estrapolazione quadratica il valore yN da yN −3 ,
yN −2 e yN −1 e lo si inserisca nell’approssimazione. Si verifichi che l’approssimazione risultante è y ′ (xN −1 ) ≈ (3yN −1 −4yN −2 +yN −3 )/(2h). Determinare
l’ordine dell’approssimazione risultante.
4. Data l’approssimazione del secondo ordine della derivata seconda y ′′ (xN −1 ) ≈
(yN − 2yN −1 + yN −2 )/h2 , si ricavi per estrapolazione cubica il valore yN da
yN −4 , yN −3 , yN −2 e yN −1 e lo si inserisca nell’approssimazione. Si verifichi
che l’approssimazione risultante è y ′′ (xN −1 ) ≈ (2yN −1 − 5yN −2 + 4yN −3 −
yN −4 )/h2 .
5. Dato il problema di trasporto

∂y
∂y


 ∂t (x, t) + c ∂x (x, t) = 0 (x, t) ∈ (a, b) × (0, +∞)
y(x, 0) = y0 (x)



y(a, t) = y0 (a)
x ∈ (a, b)
t ∈ (0, +∞)
(11.10)
87
11.3. ESERCIZI
con c > 0 e y0 (x) = y0 (a) per x ≤ a, si dimostri che la soluzione è
y(x, t) = y0 (x − ct) .
Bibliografia
[1] J.
P.
Boyd,
Chebyshev
and
Fourier
Spectral
Methods,
DOVER
Publications,
Inc.,
2000.
http://www-personal.umich.edu/~jpboyd/BOOK_Spectral2000.html
[2] , C. Canuto, M. Y. Hussaini, A. Quarteroni, and T. A. Zang, Spectral Methods
in Fluid Dynamics, Springer-Verlag, 1986.
[3] V. Comincioli, Analisi numerica: metodi, modelli, applicazioni, McGraw–Hill,
1995.
[4] E. Hairer and G. Wanner, Solving Ordinary Differential Equations II, Stiff and
Differential-Algebraic Problems, Springer, Second Revised Edition, 2002.
[5] E. Hairer, S. P. Nørsett, and G. Wanner, Solving Ordinary Differential
Equations I, Nonstiff Problems, Springer, Second Revised Edition, 2000.
[6] W.
Hundsdorfer,
Numerical
Solution
of
Advection-DiffusionReaction Equations, Lecture notes, Thomas Stieltjes Institute, 2000.
http://homepages.cwi.nl/~willem/Coll_AdvDiffReac/notes.pdf
[7] A. Iserles, A First Course in the Numerical Analysis of Differential Equations,
Cambridge Texts in Applied Mathematics, second ed., 2009.
[8] R. J. Leveque, Numerical Methods for Conservation Laws, Lectures in
Mathematics, Birkhäuser, 1992.
[9] A. Quarteroni, Modellistica numerica per problemi differenziali, Springer, terza
edizione, 2006.
88