Metodi Numerici - Equazioni Differenziali
Metodi Numerici - Equazioni Differenziali
Metodi Numerici
per le Equazioni Differenziali
a.a. 2012/13
0 Preliminari 7
5 Sistemi tridiagonali 19
6 Metodo di Newton 21
6.1 Metodo di Newton inesatto . . . . . . . . . . . . . . . . . . . 22
7 Esponenziale di matrice 23
7.1 Formula delle variazioni delle costanti . . . . . . . . . . . . . 23
7.2 Calcolo di exp(A) . . . . . . . . . . . . . . . . . . . . . . . . . 24
7.2.1 Matrici piene, di modeste dimensioni . . . . . . . . . . 24
7.2.2 Matrici sparse, di grandi dimensioni . . . . . . . . . . . 26
8 Esercizi 28
3
4 INDICE
1 BVPs 30
9 Introduzione 31
10 Differenze finite 32
10.1 Differenze finite centrate del secondo ordine . . . . . . . . . . 32
10.2 Convergenza per un problema modello . . . . . . . . . . . . . 35
10.2.1 Unicità . . . . . . . . . . . . . . . . . . . . . . . . . . . 35
10.2.2 Esistenza . . . . . . . . . . . . . . . . . . . . . . . . . 36
10.2.3 Regolarità . . . . . . . . . . . . . . . . . . . . . . . . . 36
10.2.4 Esistenza ed unicità per il problema discretizzato . . . 37
10.2.5 Consistenza . . . . . . . . . . . . . . . . . . . . . . . . 38
10.2.6 Proprietà di A . . . . . . . . . . . . . . . . . . . . . . . 38
10.2.7 Stabilità . . . . . . . . . . . . . . . . . . . . . . . . . . 38
10.2.8 Convergenza . . . . . . . . . . . . . . . . . . . . . . . . 39
10.3 Altre differenze finite . . . . . . . . . . . . . . . . . . . . . . . 40
10.3.1 Su nodi non equispaziati . . . . . . . . . . . . . . . . . 40
10.3.2 Non centrate . . . . . . . . . . . . . . . . . . . . . . . 41
10.3.3 Di ordine più elevato . . . . . . . . . . . . . . . . . . . 41
10.4 Condizioni al bordo . . . . . . . . . . . . . . . . . . . . . . . . 41
10.4.1 Condizioni di Dirichlet . . . . . . . . . . . . . . . . . . 41
10.4.2 Condizioni di Neumann . . . . . . . . . . . . . . . . . 42
10.4.3 Importanza delle condizioni al bordo . . . . . . . . . . 43
10.5 Un esempio: l’equazione della catenaria . . . . . . . . . . . . . 44
10.5.1 Iterazioni di punto fisso . . . . . . . . . . . . . . . . . . 45
10.5.2 Metodo di Newton . . . . . . . . . . . . . . . . . . . . 45
10.6 Norme ed errori . . . . . . . . . . . . . . . . . . . . . . . . . . 46
11 Metodo di shooting 48
11.1 Metodo di bisezione . . . . . . . . . . . . . . . . . . . . . . . . 48
11.2 Metodo di Newton . . . . . . . . . . . . . . . . . . . . . . . . 49
11.3 Problema ai limiti con frontiera libera . . . . . . . . . . . . . . 50
12 Equazione di Poisson 52
12.1 Equazione di Poisson bidimensionale . . . . . . . . . . . . . . 52
12.1.1 Condizioni al bordo di Dirichlet . . . . . . . . . . . . . 52
12.1.2 Condizioni al bordo miste . . . . . . . . . . . . . . . . 54
13 Metodi variazionali 56
13.1 Un problema modello . . . . . . . . . . . . . . . . . . . . . . . 56
13.1.1 Metodo di approssimazione variazionale . . . . . . . . . 58
INDICE 5
14 Esercizi 81
2 ODEs 83
15 Introduzione 84
15.1 Riduzione in forma autonoma . . . . . . . . . . . . . . . . . . 85
15.2 Equazioni di ordine superiore al primo . . . . . . . . . . . . . 85
16 Metodi ad un passo 86
16.1 Metodo di Eulero . . . . . . . . . . . . . . . . . . . . . . . . . 86
16.2 Metodo dei trapezi . . . . . . . . . . . . . . . . . . . . . . . . 88
16.3 theta-metodo . . . . . . . . . . . . . . . . . . . . . . . . . . . 90
16.3.1 Risoluzione di un metodo implicito . . . . . . . . . . . 92
16.3.2 Newton inesatto e passo variabile . . . . . . . . . . . . 94
16.3.3 Caso lineare . . . . . . . . . . . . . . . . . . . . . . . . 95
16.4 Verifica dell’implementazione . . . . . . . . . . . . . . . . . . . 95
17 Metodi multistep 97
17.1 Metodi di Adams-Bashforth . . . . . . . . . . . . . . . . . . . 97
17.2 Metodi lineari multistep . . . . . . . . . . . . . . . . . . . . . 99
17.2.1 Metodi BDF . . . . . . . . . . . . . . . . . . . . . . . . 101
17.3 Consistenza e stabilità . . . . . . . . . . . . . . . . . . . . . . 103
17.4 Influenza degli errori di arrotondamento . . . . . . . . . . . . 108
19 A-stabilità 122
19.1 A-stabilità dei metodi di Runge-Kutta espliciti . . . . . . . . . 124
19.2 A-stabilità dei metodi lineari multistep . . . . . . . . . . . . . 126
19.3 Equazioni stiff . . . . . . . . . . . . . . . . . . . . . . . . . . . 127
19.3.1 Risoluzione di un metodo implicito per problemi stiff . 129
6 INDICE
21 Esercizi 133
3 PDEs 136
22 Equazioni ADR 137
22.1 Equazione del calore . . . . . . . . . . . . . . . . . . . . . . . 137
22.1.1 Esistenza di una soluzione . . . . . . . . . . . . . . . . 137
22.1.2 Unicità della soluzione . . . . . . . . . . . . . . . . . . 140
22.2 Metodo di Fourier . . . . . . . . . . . . . . . . . . . . . . . . . 141
22.3 Metodo delle linee . . . . . . . . . . . . . . . . . . . . . . . . . 142
22.3.1 Differenze finite . . . . . . . . . . . . . . . . . . . . . . 143
22.3.2 Condizioni al bordo di Dirichlet . . . . . . . . . . . . . 144
22.3.3 Condizioni al bordo di Neumann (costanti) . . . . . . . 145
22.4 Equazione di trasporto-diffusione . . . . . . . . . . . . . . . . 145
22.4.1 Stabilizzazione mediante diffusione artificiale . . . . . . 147
22.4.2 Elementi finiti . . . . . . . . . . . . . . . . . . . . . . . 150
22.4.3 Errori spaziali e temporali . . . . . . . . . . . . . . . . 151
22.5 Esercizi . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 152
5 Bibliografia 166
Bibliografia 167
Parte 0
Preliminari
7
Capitolo 1
Interpolazione polinomiale a
tratti
8
1.1. INTERPOLAZIONE LINEARE A TRATTI 9
f (n) (ξ)
f (x) − pn−1 f (x) = (x − x1 ) · (x − x2 ) · . . . · (x − xn )
n!
ove pn−1 f è il polinomio di grado n − 1 interpolatore di f sui nodi {xi }ni=1 e ξ
un opportuno punto nell’involucro convesso di x ∪ {xi }ni=1 . Per l’interpolante
lineare a tratti pc1 f , si ha dunque per x ∈ [xi , xi+1 ]
f ′′ (ξ)
f (x) − pc1 f (x) = (x − xi )(x − xi+1 ) (1.2)
2
e pertanto
f ′′ (ξ) h2
|f (x) − pc1 f (x)| ≤ max , x ∈ [xi , xi+1 ]
ξ∈[xi ,xi+1 ] 2 4
f ′′ (ξ)
f ′ (x) − (pc1 f )′ (x) = [(x − xi+1 ) + (x − xi )]
2
e pertanto
f ′′ (ξ)
|f ′ (x) − (pc1 f )′ (x)| ≤ max 2h, x ∈ [xi , xi+1 ]
ξ∈[xi ,xi+1 ] 2
da cui
kf ′ − (pc1 f )′ k∞ ≤ hkf ′′ k∞
Capitolo 2
Formule di quadratura
gaussiana
con l’ipotesi
Z b
|x|k w(x)dx < ∞, k≥0
a
Allora, esiste un’unica famiglia {pj (x)}j , pj (x) polinomio di grado j, orto-
normale rispetto al prodotto scalare
Z b
pj (x)pi (x)w(x)dx = δij
a
ove Ln (x) è il polinomio di Lagrange che vale 1 in xn e zero in tutti gli altri
nodi, costituiscono una formula di quadratura gaussiana esatta fino al grado
polinomiale 2m − 1, cioè
Z b m
X
pj (x)w(x)dx = pj (xn )wn , 0 ≤ j ≤ 2m − 1
a n=1
10
2.1. QUADRATURA GAUSSIANA DI CHEBYSHEV(-LOBATTO) 11
In particolare
Z b m
X
δij = pj (x)pi (x)w(x)dx = pj (xn )pi (xn )wn , 0 ≤ i, j ≤ m − 1
a n=1
Nel caso in cui (a, b) sia limitato, esiste un’unica formula di quadratura esatta
fino al grado polinomiale 2m − 3 che usa come nodi x̄1 = a, x̄m = b e gli zeri
m−1
{x̄n }n=2 del polinomio di grado m − 2 della famiglia di polinomi ortogonali
rispetto alla funzione peso w(x)(x − a)(b − x). In questo caso si ha, in
particolare,
Z b m
X 0≤i≤m−3
δij = pj (x)pi (x)w(x)dx = pj (x̄n )pi (x̄n )w̄n ,
a n=1
0≤j ≤m−1
p
La famiglia {φj (x)}m
j=1 , ove φj (x) = pj−1 (x) w(x) è ovviamente ortonormale
rispetto al prodotto scalare
Z b
(f, g) = f (x)g(x)dx
a
e per essa valgono le osservazioni fatte sopra riguardo al calcolo degli integrali.
e i corrispondenti pesi
π
per n = 1 o n = m
2(m − 1)
w̄n = π
per 2 ≤ n ≤ m − 1
m−1
Capitolo 3
Ax = b (3.1)
Una strategia generale per costruire la successione {x(l) }l è basata sullo split-
ting A = P − M , ove P è non singolare. Assegnato x(1) , il termine x(l+1) è
calcolato ricorsivamente come
e(l) = Bel−1 , B = P −1 M = I − P −1 A ,
13
14 CAPITOLO 3. METODI ITERATIVI PER SISTEMI LINEARI
P z (l) = r(l)
x(l+1) = x(l) + αz (l) (3.4)
r(l+1) = r(l) − αAz (l)
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) )
3.1. METODI DI RICHARDSON 15
P z (l) = r(l)
T
z (l) r(l)
αl = T
z (l) Az (l) (3.5)
(l+1)
x = x(l) + αl z (l)
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).
Stima dell’errore
Vale la seguente stima dell’errore:
Ãp !l−1
cond (P −1 A) − 1
2
ke(l) kA ≤ 2 p ke(1) kA
−1
cond2 (P A) + 1
ke(l) k
≤ tol · cond(A)
kxk
Una modifica consiste in
che coincide con il precedente nel caso in cui come x(1) venga scelto il vettore
di zeri.
Capitolo 4
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 elemento i-esimo (i < N + 1) è la posizione
dentro data del primo elemento della colonna i e l’elemento (N +1)-esimo è il
numero totale di elementi diversi da zero incrementato di uno. 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 = [1, 3, 4, 6, 8]
17
18 CAPITOLO 4. MEMORIZZAZIONE DI MATRICI SPARSE
Sistemi tridiagonali
19
20 CAPITOLO 5. SISTEMI TRIDIAGONALI
e (
xn = yn /αn
xk = (yk − ck xk+1 )/αk , k = n − 1, n − 2, . . . 1
con un ulteriore costo O(2n) flops. GNU Octave usa automaticamente questo
algoritmo per le matrici tridiagonali.
Capitolo 6
f (x) = 0 .
kδx(r) k ≤ tol
21
22 CAPITOLO 6. METODO DI NEWTON
f = @(x) ... % f
J = @(x) ... % jacobiano di f
x0 = ... % guess iniziale
x = x0;
errest = -J(x) \ f(x);
while (norm(errest,inf) > Newt_tol)
x = x + errest;
errest = -J(x) \ f(x);
end
Esponenziale di matrice
Infatti, si ha
Z t
′ (t−t0 )a
y (t) = ae y0 +a e(t−τ )a b(τ, y(τ ))dτ +e(t−t)a b(t, y(t)) = ay(t)+b(t, y(t))
t0
23
24 CAPITOLO 7. ESPONENZIALE DI MATRICE
Si osservi che
Z t Z t
¯t
(t−τ )a 1 (t−τ )a 1 (t−τ )a ¯¯
e dτ = − −ae dτ = − e ¯ =
t0 a t0 a t0
(t−t )a
1¡ ¢ e 0
−1
= − 1 − e(t−t0 )a = (t − t0 ) =
a (t − t0 )a
= (t − t0 )ϕ1 ((t − t0 )a) ,
ove
∞
ez − 1 X z j
ϕ1 (z) = = (7.3)
z j=0
(j + 1)!
e, analogamente,
Z t
e(t−τ )a (τ − t0 )dτ = (t − t0 )2 ϕ2 ((t − t0 )a)
t0
ove
∞
ez − 1 − z X z j
ϕ2 (z) = = (7.4)
z2 j=0
(j + 2)!
A2 = (V DV −1 )2 = (V DV −1 )(V DV −1 ) = V D2 V −1
z a1 z p−1 + a2 z p−2 + . . . + ap
e ≈ , (7.5)
b1 z q−1 + b2 z q−2 + . . . + bq
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 espo-
nenziale e si fanno coincidere quanti più coefficienti possibile. Per esempio,
fissiamo p = q = 2. Si ha allora
µ ¶
z2 z3
1+z+ + + . . . (b1 z + 1) = a1 z + a2
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. Consi-
derando sempre il caso p = q = 2, si ha
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
1 ³z ´
ϕ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.
ove {di }i sono le differenze divise. Tale formula si può scrivere, nel caso
matriciale,
Ãm−1 !
Y
pm−1 (A)v = pm−2 v + dm wm , wm = (A − ξi I) v = (A − ξm−1 )wm−1
i=1
7.2. CALCOLO DI EXP(A) 27
Esercizi
28
29
BVPs
(Problemi ai limiti)
30
Capitolo 9
Introduzione
31
Capitolo 10
Differenze finite
32
10.1. DIFFERENZE FINITE CENTRATE DEL SECONDO ORDINE 33
10
errore derivata se onda
stima errore derivata se onda diff12.m
errore derivata prima
stima errore derivata prima
1
errore in norma innito
0.1
0.01
10 20 50 100
m
In Figura 10.1 si vedono gli errori (in norma infinito) tra la derivata
prima e seconda della funzione u(x) = sin(3x) e la relativa approssimazione
mediante differenze finite centrate del secondo ordine (asterischi) e le stime
h2 /6 · ku(3) k∞ e h2 /12 · ku(4) k∞ (linea continua), rispettivamente, ove h =
2π/(m − 1). In Figura 10.2 si vede invece che per la funzione u(x) = |x|7/2 ,
l’approssimazione della derivata prima mediante differenze finite centrate ha
effettivamente ordine due, mentre quella della derivata seconda no, in quanto
non esiste la derivata quarta di u(x) (h = 2/(m − 1)).
Una volta scelto il tipo di discretizzazione, invece del problema originale
(9.1) si risolve il problema discretizzato
2
∆ ui = f (xi , ui , ∆ui ),
2≤i≤m−1
u 1 = ua
u = u
m b
34 CAPITOLO 10. DIFFERENZE FINITE
1
diff12ns.m errore derivata se onda
errore derivata prima
h2
0.1
errore in norma innito
0.01
0.001
10 20 50 100
m
ove la prima e l’ultima riga devono essere trattate a parte, solitamente per
includere le condizioni al bordo. Le matrici relative alle approssimazione
della derivata prima e seconda possono essere costruite con i comandi
> toeplitz(sparse(1,2,-1/(2*h),1,m),sparse(1,2,1/(2*h),1,m));
e
> toeplitz(sparse([1,1],[1,2],[-2/h^2,1/h^2],1,m));
rispettivamente.
10.2. CONVERGENZA PER UN PROBLEMA MODELLO 35
con q, g ∈ C 0 ([a, b]), q(x) ≥ 0 per x ∈ [a, b]. La funzione q(x) dipende dal
materiale di cui è fatta la trave e g(x) è la densità di carico trasversale. La so-
luzione u(x) rappresenta il momento flettente. Vogliamo studiare l’esistenza,
l’unicità e la regolarità della soluzione analitica.
10.2.1 Unicità
Se u1 (x) e u2 (x) sono due soluzioni di (10.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 (10.3)
z(b) = 0
Poiché le funzioni integrande sono non negative, si ha che deve essere neces-
sariamente z(x) ≡ 0.
10.2.2 Esistenza
Sia z(x) = c1 z1 (x) + c2 z2 (x) la soluzione generale di −z ′′ (x) + q(x)z(x) =
0, con z1 (x) e z2 (x) indipendenti (lo spazio delle soluzioni dell’equazione
lineare omogenea ha proprio dimensione due). La soluzione di (10.3) (che
corrisponde a c1 = c2 = 0) si ottiene imponendo
(
c1 z1 (a) + c2 z2 (a) = 0
c1 z1 (b) + c2 z2 (b) = 0
cioè risolvendo un sistema lineare non singolare che ammette dunque (unica)
soluzione.
10.2.3 Regolarità
Proposizione 2. Se q, g ∈ C k ([a, b]), allora u ∈ C k+2 ([a, b]).
2 + q 2 h2 −1 0 ... ... 0
2 3
u2 g2 + ua /h2
2 3 2 3
.
..
6 7
6
6 −1 2 + q 3 h2 −1 0 ... 76 u 7 6
76 3 7 6 g3 7
7
76 . 7 6 .
.. .. .. ..
6
6 76 . 7 6 .
7
1 6 0 . . . . 0 76 . 7 6 .
7
=6
6 76 7
h2 6 .. .. .. .. 76 . 7 .
7
76 . 7 6 7 6 . 7
0 . . . . 0 76 . 7 6 .
6 7
6 7
.. 7 um−2 5 4 gm−2
6 74 5
6
4 . ... 0 −1 2 + qm−2 h2 −1 5
um−1 gm−1 + ub /h 2
0 ... ... 0 −1 2 + qm−1 h2
cioè
Au = g (10.4)
10.2.5 Consistenza
Se si sostituisce ui con la soluzione analitica u(xi ), da (10.1) si ottiene
u(xi+1 ) − 2u(xi ) + u(xi−1 ) (2)
− + q(xi )u(xi ) − g(xi ) = −τi , 2 ≤ i ≤ m − 1
h2
u(x1 ) = ua
u(xm ) = ub
10.2.6 Proprietà di A
A è 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.
10.2.7 Stabilità
Consideriamo due soluzioni relative a dati perturbati g̃ e ḡ. Si ha
Aũ = g̃
Aū = ḡ
da cui
(ũ − ū) = A−1 (g̃ − ḡ)
10.2. CONVERGENZA PER UN PROBLEMA MODELLO 39
A−1
q=0 − A
−1
= A−1
q=0 (A − Aq=0 )A
−1
≥0
kA−1 −1 T
q=0 k∞ = kAq=0 [1, . . . , 1] k∞ = max vi =
2≤i≤m−1
(b − a)2
= max v(xi ) ≤ max v(x) ≤
2≤i≤m−1 x∈[a,b] 8
10.2.8 Convergenza
Definiamo eh = [e2,h , . . . , em−1,h ]T = [u2 − u(x2 ), . . . , um−1 − u(xm−1 )]T , h =
(b − a)/(m − 1). Poiché
A[u2 , . . . , um−1 ]T = g
(2)
A[u(x2 ), . . . , u(xm−1 )]T = g − τ h
40 CAPITOLO 10. DIFFERENZE FINITE
u(x)
u3
u2
u1
u0
y = u′a (x − x2 ) + u2
h h
x0 a = x1 x2 x3
u0 = u2 − 2hu′a
in x = a sarà
u0 − 2u1 + u2 u2 − u0 u2 − 2hu′a − 2u1 + u2 u2 − (u2 − 2hu′a )
− = − =
h2 2h h2 2h
2u2 − 2u1 − 2hu′a
= − u′a = 1
h2
da cui
2u2 − 2u1 2u′a
= 1 + + u′a
h2 h
u(−1) = u(1) = 0
44 CAPITOLO 10. DIFFERENZE FINITE
0.8
errore in norma innito
u(x)
0.7 1e-4
0.6
0.5
0.4 1e-5
-1 -0.5 0 0.5 1 20 30 40 50 60
x m
10−2
fsc.m h2
errore rispetto alla soluzione esatta
errore rispetto alla soluzione di riferimento
10−3
errore in norma innito
10−4
10−5
10−6
10 102 103
m
da cui
kum − uk∞ − ε ≤ kum − um̄ k∞ ≤ kum − uk∞ + ε
se ku − um̄ k∞ = ε < kum − uk∞ . Ciò significa che si può stimare l’errore di
um usando una soluzione di riferimento um̄ solo se questa dista poco dalla
soluzione analitica e se m ≪ m̄, altrimenti la stima dice solo che kum −
um̄ k∞ / 2ε. Si ha cioè l’impressione che la soluzione numerica sia più vicina
alla soluzione analitica di quello che dovrebbe, invece è solo molto vicina a
quella di riferimento (per assurdo, se m = m̄, kum −um̄ k∞ = 0 6= kum −uk∞ ).
Una maniera molto comoda per verificare l’ordine di un metodo si basa
sulla seguente osservazione. Siano em̃ e em̂ gli errori corrispondenti a due
discretizzazioni con m̃ e m̂ punti. Supponiamo che
C
kem̃ k∞ =
(m̃ − 1)p
C
kem̂ k∞ =
(m̂ − 1)p
Si ricava
Metodo di shooting
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)
(da leggersi “valore di y1 in t dato che y2 in a vale 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 generale non lineare) F (s) = 0.
48
11.2. METODO DI NEWTON 49
∂ ∂
v(x) = u(x | u′ (a) = s) = y1 (t | y2 (a) = s)
∂s ∂s
∂ ′′ ∂
u (x) = f (x, u(x), u′ (x))
∂s ∂s
v ′′ (x) = fu (x, u(x), u′ (x))v(x) + fu′ (x, u(x), u′ (x))v ′ (x), x ∈ (a, b)
∂
v(a) = u(a | u′ (a) = s) = 0
∂s
∂ ′
v ′ (a) = u (a | u′ (a) = s) = 1
∂s
u′ (s | s + h) − u′ (s | s)
v ′ (s) = lim = −u′′ (s)
h→0 h
ove il valore u′′ (s) si ricava dal problema (11.2) e vale f (s, α, β).
Capitolo 12
Equazione di Poisson
52
12.1. EQUAZIONE DI POISSON BIDIMENSIONALE 53
uk ≈ u(xi , yj ), k = (j − 1)mx + i
u5 u6 u7 u8
u1 u2 u3 u4
ne alle differenze finite centrate del secondo ordine, senza tener conto delle
condizioni al bordo, è data da
A = Imy ⊗ Ax + Ay ⊗ Imx
ove Ax ∈ Rmx ×mx e Ay ∈ Rmy ×my . Poi, le righe di indice, diciamo k, cor-
rispondente ad un nodo al bordo vanno sostituite con il vettore della base
54 CAPITOLO 12. EQUAZIONE DI POISSON
canonica ek , diviso per h2x + h2y . Il termine noto è [b1 , b2 , . . . , bmx my ]T , ove
f (xi , yj ) se (xi , yj ) è un nodo interno, k = (j − 1)mx + i
Da (yj )
h2x +h2y se xi = a, k = (j − 1)mx + i
Db (yj )
bk = h2x +h2y se xi = b, k = (j − 1)mx + i
D (x )
h2xc+hi2y
se yj = c, k = (j − 1)mx + i
Dd (xi ) se yj = d, k = (j − 1)mx + i
h2 +h2
x y
> x = linspace(a,b,mx);
> y = linspace(c,d,my);
> [X,Y] = ndgrid(x,y);
ove
2 0 ... ... 0 2 0 ... ... 0
2 3 2 3
−2 −1
6 .7 6 .7
6 .7 6 .7
6−1
6 2 −1 0 ... .77
6−1
6 2 −1 0 ... .77
6 7 6 7
6 .. .. .. .. 7 6 .. .. .. .. 7
1 6 0 . . . . 1 6 0 . . . .
6 7 6 7
07 07
Ax = 6 7, Ay = 6 7
h2 .. .. .. .. h2 .. .. .. ..
6 7 6 7
x 6 7 y 6 7
6 0
6 . . . . 07
7
6 0
6 . . . . 07
7
6 7 6 7
6 . 7 6 . 7
6 . 7 6 . 7
4 . ... 0 −1 2 −15 4 . ... 0 −1 2 −15
0 ... ... 0 −1 2 0 ... ... 0 −2 2
(si può verificare che lo stencil [2, −5, 4, −1]/h2x è un’approssimazione al se-
condo ordine della derivata seconda). Il termine noto è [b1 , b2 , . . . , bmx my ]T ,
ove
f (xi , yj )
se (xi , yj ) è un nodo interno, k = (j − 1)mx + i
Db (yj )
h2 +h2
se xi = b, k = (j − 1)mx + i
x y
Dc (xi )
bk = 2 2 se yj = c, k = (j − 1)mx + i
hx +hy
2Na (yi )
f (xi , yj ) + hx
se xi = a, k = (j − 1)mx + i, j 6= 1, j 6= my
2Nd (xi )
f (xi , yj ) + se yj = d, k = (j − 1)my + i, i =
6 mx
hy
Capitolo 13
Metodi variazionali
ove
0
0 ≤ x < 21 − ε
1 1
gε (x) = − 2ε 2
− ε ≤ x ≤ 21 + ε
1
0 2
+ε<x≤1
56
13.1. UN PROBLEMA MODELLO 57
La “soluzione” di (13.2) è
1 1
− x 0≤x≤ −ε
2
2
µ ¶2
1 1 ε−1 1 1
uε (x) = x− + −ε≤x≤ +ε
4ε 2 4 2 2
− 1 (1 − x) 1
+ε≤x≤1
2 2
In che senso soluzione? Chiaramente u′′ε (1/2 ± ε) non esiste e quindi non è
vero che −u′′ε (x) = gε (x), x ∈ (0, 1). Cerchiamo dunque una formulazione
alternativa che renda sensato il modello per un problema cosı̀ semplice e fisico
come (13.2).
Introduciamo il seguente spazio lineare:
e il prodotto scalare su V
Z 1
(v, w) = v(x)w(x)dx
0
Qualora si trovi una soluzione debole, ha però senso verificare se per caso
non sia anche forte. Infatti, se u ∈ V è soluzione di (13.3) e u ∈ C 2 ([0, 1])
(da notare che C 2 ([0, 1]) ⊂ V ) e g è continua allora 0 = (u′ , v ′ ) − (g, v) =
(−u′′ , v) − (g, v) = −(u′′ + g, v) per ogni v ∈ V . Poiché u′′ + g è continua, si
deduce −u′′ (x) = g(x) per 0 < x < 1.
Per quanto visto, la formulazione variazionale (13.3) del problema (13.1)
è in realtà le più “fisica”: pensando al problema della trave, essa permette
di descrivere, per esempio, anche il caso in cui la densità di carico g(x)
non sia continuo. Basta infatti che sia possibile calcolare (g, v), v ∈ V e
dunque basta, per esempio, che g sia continua a tratti. Quindi, in generale,
è possibile come modello per un fenomeno fisico la sola formulazione debole.
La soluzione debole, se esiste, è unica: infatti, se u1 e u2 sono due soluzioni
di (13.3), 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
P
e l’unica possibilità per avere 0 è che wj φj (x) sia costante e dunque nullo
(poiché nullo ai bordi). Dunque, A è definita positiva.
p
ove kvk = (v ′ , v ′ ).
ku − ûk ≤ ku − vk, ∀v ∈ Vm
e quindi la tesi.
Per definizione, û è allora la proiezione ortogonale della soluzione esatta
u sul sottospazio Vm , tramite il prodotto scalare hu, vi = (u′ , v ′ ).
La scelta di Vm 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 dal fatto che se u ∈ V , allora
Z 1 Z 1 ¯1 Z 1
′ 2 ′ 2
û2 (x)dx
¯
2xû(x)û (x)dx = xû (x) dx = û (x)x¯ −
0 0 0 0
da cui p p p p
(û, û) ≤ 2 (û, û) (û′ , û′ ) = 2 (û, û) (û′ , û′ )
cioè p
(û, û) ≤ 2kûk
Poiché û soddisfa, in particolare,
da cui p
kûk ≤ 2 (g, g)
Si conclude osservando che date due perturbazioni della soluzione ũ e ū
corrispondenti rispettivamente a g̃ e ḡ, allora
e pertanto p
kũ − ūk ≤ 2 (g̃ − ḡ, g̃ − ḡ)
e cioè che piccole variazioni sui dati producono piccole variazioni sulle solu-
zioni.
φ1 φj−1 φj φj+1
h1 hj−1 hj hm−1
x1 x2 xj−2 xj−1 xj xj+1 xm
e
− 1, x <x<x
′ 1 2
φ1 (x) = h1
0, altrimenti
e
x − xm−1 , xm−1 ≤ x ≤ xm
φm (x) = hm−1
0, altrimenti
e
1
, xm−1 < x < xm
φ′m (x) = hm−1
0, altrimenti
Dunque, nell’approssimazione
m
X
û(x) = ûj φj (x)
j=1
Xm Z xi +hi
= ûj aij = g(x)φi (x)dx
j=1 xi −hi−1
Siccome il supporto di φj (x) è [xj−1 , xj+1 ], gli unici elementi non nulli aij
sono aii , ai i−1 e ai+1 i = ai i+1 . Per 1 < i < m,
Z xi µ ¶2
¶2 Z xi +hi µ
1 1 1 1
aii = (φ′i , φ′i )
= dx + − dx = +
xi −hi−1 hi−1 xi hi hi−1 hi
Z xi
1 1 1
ai i−1 = (φ′i−1 , φ′i ) = − · dx = − = ai−1 i
xi −hi−1 hi−1 hi−1 hi−1
13.1. UN PROBLEMA MODELLO 63
Per i = 1 e i = m, si ha invece
Z x1 +h1 µ ¶2
1 1
a11 = − dx =
h1 h1
Zx1x2
1 1 1
a21 = − · dx = − = a12
x2 −h1 h1 h1 h1
Z xm
1 1 1
am m−1 = − · dx = − = am−1 m
xm −hm−1 hm−1 hm−1 hm−1
Z xm µ ¶2
1 1
amm = − dx =
xm −hm−1 hm−1 hm−1
Per quanto riguarda il calcolo di (g, φi ) si può ricorrere alla formula del punto
medio: per 1 < i < m è
Z xi Z xi+1
x − xi−1 xi+1 − x
gi = (g, φi ) = g(x) dx + g(x) dx ≈
xi−1 hi−1 xi hi
µ ¶ µ ¶
xi−1 + xi hi−1 xi + xi+1 hi
≈g +g
2 2 2 2
Per i = 1 e i = m si ha invece
Z x2 µ ¶
x2 − x x1 + x2 h1
g1 = (g, φ1 ) = g(x) dx ≈ g
x1 h1 2 2
Z xm µ ¶
x − xm−1 xm−1 + xm hm−1
gm = (g, φm ) = g(x) dx ≈ g
xm−1 hm−1 2 2
L’approssimazione di
Z xi Z xi
x − xi−1
g(x)φi (x)dx = g(x) dx
xi−1 xi−1 hi−1
la formula del punto medio viene di solito sostituita dalla formula equivalente
(nel senso dell’ordine di approssimazione)
Z xi Z xi +hi
hi−1 hi
gi = (g, φi ) ≈ ḡi−1 φi (x)dx + ḡi φi (x)dx = ḡi−1 + ḡi
xi −hi−1 xi 2 2
e dunque molto simile (il termine noto è diverso, anche se dello stesso ordine)
a quella della discretizzazione con differenze finite del secondo ordine. Pertan-
to, è naturale aspettarsi, sotto opportune ipotesi di regolarità, che l’errore,
nella norma indotta dal prodotto scalare, rispetto alla soluzione analitica
tenda a zero come h2 , h = maxj hj (e ciò giustifica, a posteriori, la scelta
della formula di quadratura). Facilmente, usando la disuguaglianza (13.5) e
i risultati della Sezione 1.1.1, si può dimostrare
s
Z 1
c
ku − ûk ≤ ku − p1 uk = [u′ (x) − (pc1 u)′ (x)]2 dx ≤ hku′′ k∞
0
Da notare che il problema con due condizioni di Neumann non è ben definito,
in quanto se u(x) è soluzione, allora lo è anche u(x) + k.
Ovviamente, lo spazio Vm può essere costituito da funzioni molto più
regolari (per esempio polinomi di grado superiore).
Vediamo un approccio più generale all’implementazione del metodo degli
elementi finiti. Supponiamo di avere l elementi {ℓj }lj=1 (nel caso unidimen-
sionale, sono gli intervalli) ad ognuno dei quali sono associati due nodi. Con
ℓ1 ℓ2 ℓj−1 ℓj ℓl
x1 x2 xj−1 xj xj+1 xm
g(xℓj,1 ) + g(xℓj,2 )
gℓj = ḡj =
2
Pertanto si ha
• aij = 0, 1 ≤ i, j ≤ m, gi = 0, 1 ≤ i ≤ m
• for j = 1, . . . , l
for k = 1, 2
hj
aℓj,k ℓj,k = aℓj,k ℓj,k + φℓj,k ℓj,k , gℓj,k = gℓj,k + gℓj 2
aℓj,k ℓj,3−k = aℓj,k ℓj,3−k + φℓj,k ℓj,3−k
end
end
Se tale limite è finito per ogni k, allora la serie si dice convergere esponen-
zialmente oppure spettralmente. Significa che |uj | decade più velocemente di
13.2. METODI SPETTRALI 67
ei2π(j−k)(x−a)/(b−a)
φj (x)φk (x) =
b−a
13.2. METODI SPETTRALI 69
e quindi
Z b Z 1
ei2π(j−k)y
φj (x)φk (x)dx = (b − a)dy = 0 ,
a 0 b−a
poiché l’integrale delle funzioni sin e cos in un intervallo multiplo del loro pe-
riodo è 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
Si ha dunque
X ¡ ¢ ¡ ¢
|uj |2 = |um+1 |2 + |u0 |2 + |um+2 |2 + |u−1 |2 + |um+3 |2 + . . . =
j∈J
¡ ¢ ¡ ¢ ¡ ¢
= O (m/2)−k + O (1 + m/2)−k + O (2 + m/2)−k + . . .
ove x = (b − a)y + a.
La funzione (serie troncata di Fourier)
m m/2−1
X X
û(x) = ûj φj (x) = ûk+1+m/2 φk+1+m/2 (x) =
j=1 k=−m/2
m/2−1
X eik2π(x−a)/(b−a)
= ûk+1+m/2 √
k=−m/2
b−a
72 CAPITOLO 13. METODI VARIAZIONALI
La trasformazione
[u(x1 ), u(x2 ), . . . , u(xm )]T → [û1 , û2 , . . . , ûm ]T
si chiama trasformata di Fourier discreta √ di u e û1 , . . . , ûm coefficienti di Fou-
T
rier di u. Il vettore m · [û ,
√1 2 û , . . . , û m ] / b − a può essere scritto come pro-
imπy1
dotto matrice-vettore F m[u(x1 )e , u(x2 )eimπy2 , . . . , u(xm )eimπym ]T , ove
µ ¶r
e−i(j−1)2πyk xk − a b−a
F = (fjk ), fjk = √ = φj+m/2 .
m b−a m
Alternativamente, si può usare la Fast Fourier Transform (FFT). Il comando
fft applicato al vettore [u(x1 )eimπy √ , u(x2 )e
1 imπy2
, . . . , u(xm )eimπym ]T produ-
T
ce il vettore m · [û1 , û2 , . . . , ûm ] / b − a, cosı̀ come il comando fftshift
applicato al risultato del comando fft applicato a [u(x1 ), u(x2 ), . . . , u(xm )].
Dati dei coefficienti v̂j , j = 1, . . . , m, si può considerare la funzione
(periodica)
X m
v̂j φj (x)
j=1
La trasformazione
√
Il vettore b − a · [v̂ˆ1 eimπy1 , v̂ˆ2 eimπy2 , . . . , v̂ˆm eimπy√
m T
] /m può essere scritto
′ T
come prodotto matrice-vettore F [v̂1 , v̂2 , . . . , v̂m ] / m, ove F ′ denota, co-
me in GNU Octave, la trasposta coniugata di F . Alternativamente, √ il co-
mando ifft applicato al vettore [v̂1 , v̂2 , . . . , v̂m ] produce il vettore b − a ·
[v̂ˆ1 eimπy1 , v̂ˆ2 eimπy2 , . . . , v̂ˆm eimπym ]/m, mentre, se applicato al risultato del co-
mando
√ ifftshift applicato al vettore [v̂1 , v̂2 , . . . , v̂m ], produce il vettore
b − a · [v̂ˆ1 , v̂ˆ2 , . . . , v̂ˆm ]/m.
Da notare che F ′ = F −1 .
10−2
ordine 1
ordine 2
fourier.m
ordine 4
ordine 8
10−4 errore
10−6
10−8
errore in norma-2
10−10
10−12
10−14
10−16
10−18
1 10 100
ei(j−1−m/2)2π(x−a)/(b−a)
φj (x) = √
b−a
da cui
m
X Z π m
X Z π Z π
− uj φ′′j (x)φk (x)dx + uj φj (x)φk (x)dx = g(x)φk (x)dx
j=1 −π j=1 −π −π
Poiché µ ¶2
i(j − 1 − m/2)2π
φ′′j (x) = φj (x) = λ2j φj (x)
b−a
usando l’ortonormalità delle funzioni di Fourier e calcolando i coefficienti di
Fourier di g(x), si ha
da cui
ĝk
ûk = , 1≤k≤m
1 − λ2k
e quindi
m
X
u(x) ≈ ûj φj (x)
j=1
1
http://math-atlas.sourceforge.net/
2
http://xianyi.github.com/OpenBLAS/
76 CAPITOLO 13. METODI VARIAZIONALI
Oppure si può costruire la matrice F relativa ai nodi (ciò funziona anche per
nodi non equispaziati). Infine, si può usare la trasformata di Fourier non
equispaziata NFFT.
ove gli {xn } sono i nodi della quadratura gaussiana relativa alla famiglia {φj }
è la stessa del problema di Galerkin
m
X Z b Z b
ûj Lφj (x)φi (x)w(x)dx = g(x)φi (x)w(x)dx
j=1 a a
m
à m ! m
à m !
X X X X
ûj Lφj (xn )φi (xn )wn = ûj Lφj (xn )φi (xn )wn =
n=1 j=1 j=1 n=1
Xm
= g(xn )φi (xn )wn , 1≤i≤m
n=1
Anche in questo caso il metodo di collocazione può essere visto come un me-
todo di Galerkin pseudospettrale: basta prendere come nodi di collocazione
gli m − 2 nodi di quadratura gaussiana. Si ha poi
m Ãm−2 ! m−2
X X X
û Lφ (x )φ (x )w = g(xn )φi (xn )wn , 1 ≤ i ≤ m − 2
j j n i n n
j=1
n=1 n=1
X m
ûj φj (a) = ua
j=1
m
X
ûj φ′j (b) = u′b
j=1
Collocazione Gauss–Lobatto–Chebyshev
I polinomi di Chebyshev sono definiti da
Tj (x) = cos(j arccos(x)), −1 ≤ x ≤ 1
e soddisfano
ππ i=j=0
Z
1
Tj (x)Ti (x)
√ dx = i = j 6= 0
−1 1 − x2
2
0 i 6= j
(lo si vede con il cambio di variabile x = cos θ e applicando le formule di Wer-
ner, oppure integrando per parti due volte). I nodi di (Gauss–)Chebyshev–
Lobatto sono xn = cos((n − 1)π/(m − 1)), 1 ≤ n ≤ m. Possiamo allora
definire la seguente famiglia di funzioni ortonormali
r r
1 2
φ1 (x) = T0 (x), φj (x) = Tj−1 (x), 2 ≤ j ≤ m
π π
13.3. METODI DI COLLOCAZIONE 79
bordo) è
q(x1 ) u1 g(x1 )
T
q(x2 ) u2 g(x2 )
T
−T2 + ... T0 .. = ..
. .
q(xm ) um g(xm )
Esercizi
81
82 CAPITOLO 14. ESERCIZI
u(0) = 0
u(1) = 0
8. Noti gli zeri dei polinomi di Legendre e i pesi di quadratura della ri-
spettiva formula gaussiana, si ricavino i nodi e i pesi di una formula
gaussiana nell’intervallo [a, b] rispetto al peso w(x) = 1.
ODEs
(Equazioni differenziali
ordinarie)
83
Capitolo 15
Introduzione
84
15.1. RIDUZIONE IN FORMA AUTONOMA 85
Metodi ad un passo
86
16.1. METODO DI EULERO 87
si vede che l’errore locale coincide con y(tn+1 ) − y ∗n+1 , cioè con la differenza
tra la soluzione esatta al tempo tn+1 e la soluzione approssimata al tempo
tn+1 che si otterrebbe applicando il metodo numerico al problema differenziale
e supponendo esatta la soluzione al tempo tn . È chiaro allora che ad ogni
passo si commette un nuovo errore che si accumula con l’errore prodotto ai
passi precedenti. Per il metodo di Eulero, si ha
ove en,k = y n,k − y(tn ). La quantità maxn ken,k k si chiama errore globale.
si ricava
en+1,k = en,k + k[f (tn , y n ) − f (tn , y(tn ))] + O(k 2 )
da cui, siccome f è lipschitziana,
Allora
c
k[(1 + kλ)n − 1], 0 ≤ n ≤ m
ken,k k ≤
λ
(si dimostra per induzione). Poiché 1 + kλ < ekλ , (1 + kλ)n < enkλ ≤ emkλ =
∗
et λ . Dunque
c ∗
ken,k k ≤ k(et λ − 1), 0 ≤ n ≤ m
λ
da cui
lim+ max ken,k k = 0
k→0 0≤n≤m
k
y(tn+1 ) − y(tn ) − (f (tn , y(tn )) + f (tn+1 , y(tn+1 ))) =
2
2
k k3 k
= y(tn )+ky ′ (tn )+ y ′′ (tn )+ y ′′′ (tn )+O(k 4 )−y(tn )− (y ′ (tn )+y ′ (tn+1 )) =
2 6 2
2 3
k k
ky ′ (tn ) + y ′′ (tn ) + y ′′′ (tn ) + O(k 4 )+
µ 2 6 ¶
k k 2 ′′′
− ′ ′ ′′
y (tn ) + y (tn ) + ky (tn ) + y (tn ) + O(k ) = O(k 3 )
3
2 2
(supponendo y ∈ C 3 ). Dunque il metodo dei trapezi è di ordine 2.
Teorema 7. Il metodo dei trapezi è convergente.
Dimostrazione. Si ha
k
en+1,k = en,k + (f (tn , y n ) − f (tn , y(tn )))+
2
k
+ (f (tn+1 , y n+1 ) − f (tn+1 , y(tn+1 ))) + O(k 3 )
2
da cui
kλ
ken+1,k k ≤ ken,k k + (ken,k k + ken+1,k k) + ck 3 , c>0
2
Se k < 2/λ,
à ! à !
kλ
1+ c
ken+1,k k ≤ 2
kλ
ken,k k + k3
1− 2
1 − kλ
2
Allora "Ã !n #
kλ
ck 2 1+ 2
ken,k k ≤ kλ
−1 , 0≤n≤m
λ 1− 2
90 CAPITOLO 16. METODI AD UN PASSO
100
pendoloperiod.m
10−1
m
10−2
4 periodi
error in innity norm
3 periodi
10−3 2 periodi
10−4 1 periodo
m−2
10−5
10−6
200 400 600 800
m
Figura 16.1: Errori del metodo dei trapezi per l’equazione del pendolo
linearizzata.
In Figura 16.1 si vede chiaramente che l’errore (per l’equazione del pen-
dolo linearizzata) decade come k 2 ∝ m−2 e che lo stesso cresce linearmente
con l’intervallo di tempo. Entrambi i metodi descritti sono ad un passo (cioè
la soluzione y n+1 dipende esplicitamente solo da y n ). Il metodo dei trapezi è
però implicito, cioè la soluzione y n+1 è implicitamente definita dall’equazione
(in generale non lineare)
k k
Fn (y n+1 ) = y n+1 − f (tn+1 , y n+1 ) − y n − f (tn , y n ) = 0
2 2
16.3 θ-metodo
Il θ-metodo è una generalizzazione dei metodi precedenti e si scrive
y n+1 = y n + k[(1 − θ)f (tn , y n ) + θf (tn+1 , y n+1 )], n≥0
(16.4)
y 0 = y(t0 )
16.3. THETA-METODO 91
da cui
y(tn+1 ) − y ∗n+1 = kθ[f (tn+1 , y(tn+1 )) − f (tn+1 , y ∗n+1 )]+
µ ¶ µ ¶
1 2 ′′ 1 θ
+ − θ k y (tn ) + − k 3 y ′′′ (tn ) + O(k 4 )
2 6 2
La funzione g soddisfa
ed è una contrazione se
1
kθλ < 1 ⇒ k < (16.6)
θλ
Dunque, a patto di prendere k sufficientemente piccolo, si può ricavare la
soluzione x = y n+1 con il metodo del punto fisso, partendo, per esempio,
da x(0) = y n . Il metodo del punto fisso, facile da applicare, ha convergenza
lineare ed è inoltre soggetto ad una restrizione sul passo k che potrebbe essere
eccessiva. Per questi motivi è utile considerare anche il metodo di Newton
per la risoluzione del sistema non lineare. La matrice jacobiana associata è
µ ¶
∂fi (tn+1 , x)
Jn (x) = I − kθ
∂xj ij
• r=0
• x(r) = y n
end
• y n+1 = x(r)
La tolleranza Newt tol va presa tenendo conto che si sta comunque commet-
tendo un errore proporzionale a k 2 (trapezi) o addirittura k. Come valore
iniziale x(0) , invece di y n si può prendere y n + kf (tn , y n ), che corrisponde
a y(tn+1 ) approssimato con il metodo di Eulero e dunque dovrebbe essere
un’approssimazione migliore di y n+1 (ma si veda il capitolo 19). Anche altre
scelte possono andare bene, a patto di avere convergenza del metodo.
94 CAPITOLO 16. METODI AD UN PASSO
• x(r) = y n
• Pn Jn (x(r) ) = Ln Un
end
• y n+1 = x(r)
log em2 ,k2 − log em1 ,k1 = p(log k2 − log k1 ) = −p(log m2 − log m1 ).
Metodi multistep
97
98 CAPITOLO 17. METODI MULTISTEP
Dunque possono essere calcolati una volta per tutte. Calcoliamo l’ordine di
tale metodo: come al solito, dobbiamo valutare l’espressione
s−1
X
y(tn+s ) − y(tn+s−1 ) − k bj f (tn+j , y(tn+j ))
j=0
ky (s+1) (τ̄ )k
kf (τ, y(τ )) − q(τ )k ≤ |(τ − tn ) · . . . · (τ − tn+s−1 )| ≤
s! | {z }
s termini
(s+1)
ky (τ̄ )k
≤ s!k s = O(k s ), tn+s−1 ≤ τ ≤ tn+s
s!
e quindi
s−1
X
y(tn+s ) − y(tn+s−1 ) − k bj f (tn+j , y(tn+j )) =
j=0
Z tn+s s−1
X
f (τ, y(τ ))dτ − k bj f (tn+j , y(tn+j )) =
tn+s−1 j=0
Z tn+s Z tn+s
f (τ, y(τ ))dτ − q(τ )dτ = O(k s+1 )
tn+s−1 tn+s−1
y n+1 = y n + kf (tn , y n )
17.2. METODI LINEARI MULTISTEP 99
da cui
k 3k
y n+2 = y n+1 − f (tn , y n ) + f (tn+1 , y n+1 ) (17.3)
2 2
Il valore y 1 può essere ricavato, per esempio, anche con il metodo di Eulero,
poiché si ha, in tal caso, y 1 = y(t1 ) + O(k 2 ).
y n+2 − 3y n+1 + 2y n =
k £ ¤
= −5f (tn , y n ) − 20f (tn+1 , y n+1 ) + 13f (tn+2 , y n+2 ) (17.7)
12
si ha
a0 + a1 + a2 =0
a1 + 2a2 = b0 + b1 + b2 ⇒ ordine (almeno) 1
a1 + 4a2 = 2(b1 + 2b2 ) ⇒ ordine (almeno) 2
a1 + 8a2 6= 3(b1 + 4b2 ) ⇒ ordine 2
e kbs f (tn+s , y n+s ) ≈ kbs y ′ (tn + sk), conviene cercare una combinazione li-
neare di y n+j , 0 ≤ j ≤ s che approssimi kbs y ′ (tn+s ). Si procede dunque con
102 CAPITOLO 17. METODI MULTISTEP
(
y(tn ) = y(tn + k) − ky ′ (tn + k) + O(k 2 )
y(tn + k) = y(tn + k)
4k 2 ′′
y(tn ) = y(tn + 2k) − 2ky ′ (tn + 2k) + y (tn + 2k) + O(k 3 )
2
′ k 2 ′′
y(tn + k) = y(tn + 2k) − ky (tn + 2k) + y (tn + 2k) + O(k 3 )
2
y(tn + 2k) = y(tn + 2k)
da cui
1
a0 =
a 2 = 1
3
a = −4
a0 + a1 + a2 = 0
1
− 2a − a = b ⇒ 3 (17.8)
0 1 2
a = 1
2a0 + a1 = 0
2
2
2
b2 =
3
e il metodo è di ordine 2.
I metodi BDF sono gli unici metodi multistep in cui non è difficile calcolare
i coefficienti anche nel caso di passi temporali variabili. Sempre per s = 2,
se tn+1 = tn + kn+1 e tn+2 = tn+1 + kn+2 , allora
y(tn ) = y(tn + kn+1 + kn+2 ) − (kn+1 + kn+2 )y ′ (tn + kn+1 + kn+2 )+
(kn+1 + kn+2 )2 ′′
+ y (tn + kn+1 + kn+2 ) + . . .
2
y(tn + kn+1 ) = y(tn + kn+1 + kn+2 ) − kn+2 y ′ (tn + kn+1 + kn+2 )+
2
kn+2
+ y ′′ (tn + kn+1 + kn+2 ) + . . .
2
y(tn + kn+1 + kn+2 ) = y(tn + kn+1 + kn+2 )
17.3. CONSISTENZA E STABILITÀ 103
da cui i coefficienti
2
kn+2
a 0 =
a2 = 1
(kn+1 + 2kn+2 )kn+1
2
a = − (kn+1 + kn+2 )
a0 + a1 + a2 = 0
1
− a0 (kn+1 + kn+2 ) − a1 kn+2 = b2 kn+2 ⇒ (kn+1 + 2kn+2 )kn+1
0 (kn+1 + kn+2 )2 + a1 k 2 = 0
a
a2 = 1
2 2 n+2
b2 =
(kn+1 + kn+2 )
(kn+1 + 2kn+2 )
yn = e + (2n−1 − 1)ε
104 CAPITOLO 17. METODI MULTISTEP
come ben si vede dalla Tabella 17.1. Abbiamo quindi un metodo la cui
soluzione numerica diverge facendo tendere il passo temporale a 0 (cioè pro-
prio l’opposto di quanto dovrebbe succedere). È proprio un piccolo errore
commesso ad un passo che si accumula in maniera distruttiva. Si potrebbe
pensare che il problema nasca dall’arrotondamento dell’aritmentica di mac-
china e che in assenza di questo tutto possa funzionare, ma non è cosı̀. Per
il problema (17.9) consideriamo un metodo ad s passi di ordine p i cui primi
valori siano
zj = k p θj + e, 0 ≤ j ≤ s − 1
ove θ è una radice del polinomio ρ(w) (che d’ora in poi chiameremo carat-
teristico). Sono valori accettabili, in quanto distano O(k p ) dalla soluzione
analitica (per contro, e + ε dista O(1) per qualunque k). Ora zs si trova
risolvendo s
X
aj zj = 0
j=0
da cui
s−1
X s−1
X s−1
X
p j
zs = − aj zj = −k aj θ − e aj = k p θ s + e
j=0 j=0 j=0
P
ove si è usato ρ(θ) = 0 e j aj = 0. Per induzione si arriva a provare che
zn,k = zn = k p θn + e
17.3. CONSISTENZA E STABILITÀ 105
|zn,k − yn,k | = k p θn
mentre
|zn,k − yn,k | = k p nθn
e quindi le soluzioni divergono se |θ| ≥ 1. È giustificata allora (in analogia
con quanto visto al paragrafo 10.2.7) la seguente
106 CAPITOLO 17. METODI MULTISTEP
z ij = y j + δ ij , 0≤j ≤s−1
X s s
X
aj z in+j = k bj f (tn+j , z in+j ) + δ in+s , 0 ≤ n ≤ m − s
j=0 j=0
Dunque, perché un metodo sia convergente (cioè l’errore tenda a zero con
k), occorre che le perturbazioni della soluzione introdotte ad ogni passo (da
errori di approssimazione del metodo stesso o arrotondamento) rimangano
limitate (ciò non fa esplodere la prima parte di (17.10)) e che l’errore locale
tenda a zero con k (ciò fa andare a zero la seconda parte di (17.10)).
Abbiamo visto che la condizione delle radici è necessaria affinché pertur-
bazioni limitate generino soluzioni limitate. In realtà essa è anche sufficiente.
Inoltre questa condizione permette non solo di mantenere limitate le pertur-
bazioni, ma anche di farle tendere a 0. Si ha infatti il seguente teorema
fondamentale:
Teorema 10 (Equivalenza di Dahlquist). Un metodo lineare multistep con
valori iniziali consistenti è convergente se e solo se è consistente e stabile
(cioè il suo polinomio caratteristico soddisfa la condizione delle radici).
La grande portata di questo teorema è che il risultato è valido non solo
per il problema modello (17.9). Inoltre, se i valori iniziali approssimano con
ordine p le soluzioni analitiche, allora il metodo è convergente con ordine
p. Ritornando al metodo (17.7), si ha che θ = 2 è radice del polinomio
1010
k = 0.125
k = 0.250
msinstabile.m
k = 0.500
k = 0.625
soluzione esatta
108
106
104
102
0 1 2 3 4 5
t
che abbiamo usato come modello. Pertanto, il metodo non è stabile, come
si vede anche in Figura 17.2 ove il metodo è stato applicato al problema
differenziale (
y ′ (t) = y(t), t ∈(0, 5]
y(0) = 1
Come corollario al teorema precedente, abbiamo che ogni metodo ad un
passo è stabile (perché ρ(w) = w − 1) e che i metodi di Adams–Bashforth
e Adams–Moulton sono stabili (perché ρ(w) = ws − ws−1 ). Esiste un limite
superiore per l’ordine di un metodo a s passi, dato dal seguente
ck 2 + ε
ken k ≤ [(1 + kλ)n − 1] + (1 + kλ)n ke0 k, 0≤n≤m
kλ
ove ε = max0≤n≤m εn e quindi
µ ¶
t∗ λ ck ε ∗
ken k ≤ (e − 1) + + et λ ke0 k, 0≤n≤m
λ kλ
17.4. INFLUENZA DEGLI ERRORI DI ARROTONDAMENTO 109
Metodi di Runge–Kutta
110
18.1. METODI DI RUNGE–KUTTA ESPLICITI 111
cioè
ξ 1 = y n ≈ y(tn )
µ ¶
k k
ξ 2 = y n + f (tn , ξ 1 ) ≈ y tn +
2 2
µ ¶
k
y n+1 = y n + kf tn + , ξ 2 ≈ y(tn+1 )
2
L’idea generale dei metodi espliciti di Runge–Kutta è quella, come al solito,
di sostituire l’integrale nella formula risolutiva
Z tn+1
y(tn+1 ) = y(tn ) + f (τ, y(τ ))dτ
tn
ν
X
y(tn+1 ) ≈ y(tn ) + k bj f (tn + cj k, y(tn + cj k))
j=1
0 (0)
c2 a2,1
c3 a3,1 a3,2
.. .. .. ...
. . .
cν−1 aν−1,1 aν−1,2 ... aν−1,ν−2
cν aν,1 aν,2 ... aν,ν−2 aν,ν−1
b1 b2 ... bν−2 bν−1 bν
Tabella 18.1: Tableau di Butcher per metodi di Runge–Kutta espliciti.
y(tn+1 )−y(tn )−kb1 f (tn , y(tn ))−kb2 f (tn +c2 k, y(tn )+a2,1 kf (tn , y(tn ))) =
k2
= y(tn ) + ky ′ (tn ) + y ′′ (tn ) + O(k 3 ) − y(tn ) − kb1 y ′ (tn )+
· 2 ¸
∂f ∂f ′ 2
−kb2 f (tn , y(tn )) + (tn , y(tn ))c2 k + (tn , y(tn ))a2,1 ky (tn ) + O(k ) =
∂t ∂y
k2
= ky ′ (tn ) + y ′′ (tn ) + O(k 3 ) − kb1 y ′ (tn ) − kb2 y ′ (tn )+
2 · ¸
∂f 2 ∂f ′
− kb2 (tn , y(tn ))c2 k − k a2,1 b2 (tn , y(tn ))y (tn ) =
∂t ∂y
k 2 ′′
′
= ky (tn ) + y (tn ) + O(k 3 ) − kb1 y ′ (tn ) − kb2 y ′ (tn )+
2 · ¸
∂f 2 ′′ ∂f
− kb2 (tn , y(tn ))c2 k − k a2,1 b2 y (tn ) − (tn , y(tn )) =
∂t ∂t
µ ¶
′ 2 1
= k(1 − b1 − b2 )y (tn ) + k − a2,1 b2 y ′′ (tn )+
2
∂f
− k 2 (b2 c2 − a2,1 b2 ) (tn , y(tn )) + O(k 3 )
∂t
0 0 0
1 1 2 2
1 1 2 2 3 3
1 1 1 3
2 2
0 1 4 4
numero stadi ν 1 2 3 4 5 6 7 8
massimo ordine p 1 2 3 4 4 5 6 6
c1 a1,1
c2 a2,1 a2,2
c3 a3,1 a3,2 a3,3
.. .. .. ... ...
. . .
cν−1 aν−1,1 aν−1,2 ... aν−1,ν−2 aν−1,ν−1
cν aν,1 aν,2 ... aν,ν−2 aν,ν−1 aν,ν
b1 b2 ... bν−2 bν−1 bν
ν
X
ai,j = ci , 1≤i≤ν
j=1
e
ν
X
bj = 1
j=1
18.2. UN ESEMPIO DI METODO DI RUNGE–KUTTA IMPLICITO 115
da cui si deduce che (y n +y n+1 )/2 = ξ 1 e dunque (18.6b) coincide con (18.7).
È un metodo di ordine 2 (lo si dimostra a partire da (18.7) con passaggi molto
simili a quelli fatti per i metodi di Runge–Kutta espliciti di ordine 2) e in
qualche modo simile al metodo dei trapezi. Gode della seguente importante
proprietà: se per
y ′ (t) = f (t, y(t))
y(t)T y(t) è costante, cioè xT f (t, x) = 0 per ogni t, allora per (18.7) vale
yT T
n+1 y n+1 = y n y n . Prima di dimostrarlo, osserviamo che la proprietà è inte-
ressante quando y(t) è un vettore di dimensione maggiore di uno, altrimenti
deve essere f = 0 e dunque banalmente y n+1 = y n . Osserviamo poi che
µ ¶
k y + y n+1
y n+1 − y n = kf tn + , n
2 2
da cui la tesi.
3.3
symplectic.m Heun
trapezi
punto medio impli ito
3.25
3.2
3.15
y1 (t)2 + y2 (t)2 + y3 (t)2
3.1
3.05
2.95
2.9
2.85
2.8
0 20 40 60 80 100 120
libero):
′ I2 − I3
y1 (t) =
y2 (t)y3 (t)
I2 I3
I3 − I1
y2′ (t) = y3 (t)y1 (t) (18.8)
I3 I1
I1 − I2
y3′ (t) =
y1 (t)y2 (t)
I1 I2
y rappresenta il momento angolare e I1 , I2 e I3 sono i momenti principali
d’inerzia. Evidentemente |y(t)|2 = y1 (t)2 + y2 (t)2 + y3 (t)2 è costante. Pro-
viamo a risolvere il sistema con i metodi punto medio implicito, trapezi e
Heun (tutti di ordine 2): l’integrazione fino al tempo t∗ = 120 con 300 passi
temporali (I1 = 2, I2 = 1, I3 = 2/3, y(0) = [1, 1, 1]T ) produce il grafico
in Figura 18.1, in cui si vede che il metodo punto medio implicito conserva
“esattamente” la quantità |y(t)|2 , a differenza dei metodi dei trapezi e di
Heun. Ovviamente bisogna tenere in conto che gli errori di approssimazione
(risoluzione del sistema non lineare) non garantiscono l’esatta uguaglianza
yT T
n+1 y n+1 = y n y n anche per il metodo punto medio implicito. I tre metodi
testati sono dello stesso ordine, ma uno produce soluzioni “qualitativamente”
migliori.
È buona norma, dopo aver definito la matrice A e i vettori c e b del
tableau, controllare che queste condizioni siano soddisfatte, confrontando
sum(A,2) con c e sum(b) con 1, al fine di evitare banali errori. Ovviamente,
si tratta di condizioni necessarie ma non sufficienti per garantire la corretta
implementazione.
Anche il θ-metodo può essere fatto rientrare nella classe dei metodi di
Runge–Kutta semiimpliciti:
ξ1 = yn
ξ 2 = y n + k(1 − θ)f (tn , ξ 1 ) + kθf (tn + k, ξ 2 )
y
n+1 = y n + k(1 − θ)f (tn , ξ 1 ) + kθf (tn + k, ξ 2 )
0
c2 a2,1
c3 a3,1 a3,2
.. .. .. ...
. . .
cν−1 aν−1,1 aν−1,2 ... aν−1,ν−2
b1 b2 ... bν−2 bν−1
0
c2 a2,1
c3 a3,1 a3,2
.. .. .. ...
. . .
cν−1 aν−1,1 aν−1,2 ... aν−1,ν−2
cν aν,1 aν,2 ... aν,ν−2 aν,ν−1
b̂1 b̂2 ... b̂ν−2 b̂ν−1 b̂ν
Tabella 18.5: Metodi di Runge–Kutta di ordine p − 1 e p.
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
dice embedded e si scrive di solito un unico tableau, come nella Tabella 18.6.
Il fatto che per trovare metodi di Runge–Kutta sia necessario risolvere si-
stemi non lineari per i coefficienti, rende difficile ma non impossibile trovare
coppie di metodi con tali caratteristiche.
Consideriamo il sistema differenziale
( ′
ỹ (t) = f (t, ỹ(t))
ỹ(tn ) = y (p)
n
(p)
ove y n è l’approssimazione di y(tn ) ottenuta con il metodo di Runge–Kutta
di ordine p. Si ha allora
(p) (p−1) (p) p (p−1)
ky n+1 − y n+1 k = ky n+1 − ỹ(tn+1 ) + ỹ(tn+1 ) − y n+1 k ≤ Cn+1 kn+1 , (18.9)
18.3. METODI DI RUNGE-KUTTA EMBEDDED 119
0
c2 a2,1
c3 a3,1 a32
.. .. .. ...
. . .
cν aν,1 aν,2 ... aν,ν−1
b1 b2 ... bν−1
b̂1 b̂2 ... b̂ν−1 b̂ν
Tabella 18.6: Metodi di Runge–Kutta embedded di ordine p − 1 e p.
p
(cioè C̃n+1 = Cn+1 ) e si impone che l’errore C̃n+1 k̃n+1 valga proprio quanto
la tolleranza richiesta, ricavando
à (p−1)
!1/p à (p−1)
!1/p
tola + ky n+1 k · tolr tola + ky n+1 k · tolr
k̃n+1 = = (p) (p−1)
· kn+1
C̃n+1 ky n+1 − y n+1 k
(1) (2)
• y n+1 = y n + kn+1 f 1 (metodo di Eulero)
(2−1)
• en+1 = kn+1 [(b̂1 − 1)f 1 + b̂2 f 2 ]
(2−1) (1)
• IF ken+1 k > tola + ky n+1 ktolr
n = n − 1 (time step rifiutato)
ELSE
(2) (1) (2−1)
y n+1 = y n+1 + en+1
END
h i1/2
(1) (2−1)
• kn+2 = (tola + ky n+1 k · tolr )/ken+1 k · kn+1
• n=n+1
Forse il più importante metodo di Runge–Kutta embedded è il Runge–
Kutta–Fehlberg, di ordine (4)5, il cui tableau è riportato in Tabella 18.7.
18.3. METODI DI RUNGE-KUTTA EMBEDDED 121
0
1 1
4 4
3 3 9
8 32 32
12 1932
13 2197
− 7200
2197
7296
2197
439 3680 845
1 216
−8 513
− 4104
1 8 3544 1859
2
− 27 2 − 2565 4104
− 11
40
25 1408 2197
216
0 2565 4104
− 15
16 6656 28561 9 2
135
0 12825 56430
− 50 55
A-stabilità
la cui soluzione esatta y(t) = eλ(t−t0 ) y0 tende a zero, per t → +∞, se ℜ(λ) <
0. Analizziamo il comportamento del metodo di Eulero per questo problema,
supponendo di avere fissato il passo temporale k: si ha
da cui
yn = (1 + kλ)n y0
Si ha
da cui
2ℜ(λ)
lim yn = 0 ⇔ k < − (19.2)
n→∞ |λ|2
Dunque, la soluzione numerica ottenuta con il metodo di Eulero ha lo stesso
comportamento della soluzione analitica solo se il passo temporale è sufficien-
temente piccolo. Altrimenti, la soluzione può essere completamente diversa
(limn→∞ |yn | = |y0 | o limn→∞ yn = ∞). Nel caso di Eulero implicito, invece,
si ha µ ¶n
1
yn = y0
1 − kλ
122
123
da cui
disuguaglianza sempre soddisfatta, poiché ℜ(λ) < 0. Anche per il metodo dei
trapezi la soluzione numerica tende a 0 per n → ∞. Ma non è vero, in gene-
rale, per qualunque metodo implicito. Analizziamo infatti il comportamento
generale del θ-metodo per questo problema: si ha
da cui · ¸n
1 + (1 − θ)kλ
yn = y0
1 − θkλ
Si ha
¯ ¯
¯ 1 + (1 − θ)kλ ¯
lim yn = 0 ⇔ ¯¯ ¯ < 1 ⇔ |1 + (1 − θ)kλ| < |1 − θkλ| ⇔
n→∞ 1 − θkλ ¯
0 < (θ2 − (1 − θ)2 )k 2 ℜ(λ)2 − (2θ + 2(1 − θ))kℜ(λ) + (θ2 − (1 − θ)2 )k 2 ℑ(λ)2
da cui
lim yn = 0 ⇔ 0 < (2θ − 1)k 2 |λ|2 − 2kℜ(λ)
n→∞
2ℜ(λ)
lim yn = 0 ⇔ k < , (2θ − 1 < 0) (19.3)
n→∞ (2θ − 1)|λ|2
4 4 4
2 2 2
0 0 0
-2 -2 -2
-4 -4 -4
-4 -2 0 2 4 -4 -2 0 2 4 -4 -2 0 2 4
si ha ¯ ¯
¯θ − 1¯
lim |r(kλ)| = ¯¯ ¯
kℜ(λ)→−∞ θ ¯
Tale limite vale proprio 1 per θ = 1/2. Significa che se ℜ(λ) ≪ 0 oppure k ≫
0 il metodo dei trapezi potrebbe mostrare qualche problema di instabilità. In
tal caso, il metodo migliore, da questo punto di vista, è il metodo di Eulero
implicito (θ = 1). In Figura 19.2 vediamo l’applicazione dei due metodi al
problema (
y ′ (t) = −2000(y − cos t), t ≤ 1.5
(19.4)
y(0) = 0
Se
lim |r(kλ)| = 0
kℜ(λ)→−∞
2
Trapezi
Eulero impli ito
Lstability.m
1.5
1
y(t)
0.5
-0.5
0 0.5 1 1.5 2
Figura 19.2: Metodi dei trapezi e di Eulero implicito per la soluzione di (19.4)
con k = 1.5/40.
4 4 4
2 2 2
0 0 0
-2 -2 -2
-4 -4 -4
-4 -2 0 2 4 -4 -2 0 2 4 -4 -2 0 2 4
da cui la tesi.
Dimostrazione. Si ha
la condizione (19.3) per il metodo di Eulero impone k < 1/50 = 0.02. D’altra
parte, la soluzione analitica del problema per t∗ = 0.4 è minore di 10−17 (e
dunque, trascurabile, nel senso che y(0) − t(t∗ ) = y(0), in precisione doppia).
Dunque, con poco più di 20 passi il metodo di Eulero arriva a calcolare
adeguatamente la soluzione sino a t∗ .
10300
Eulero espli ito
Eulero impli ito stiff.m
10250
10200
errore in norma innito
10150
10100
1050
100
k < 2/100
10−50
500 1000 2000 2500
m
Figura 19.4: Eulero esplicito e Eulero implicito per la soluzione di (19.5) fino
al tempo t∗ = 40.
e la sua norma infinito è minore di 10−17 per t∗ = 40. Poiché però per poter
calcolare la prima componente serve un passo temporale k < 0.02, sono ne-
cessari più di 2000 passi (vedi Figura 19.4), anche se la prima componente
diventa trascurabile dopo pochi passi e la seconda non richiederebbe un cosı̀
elevato numero di passi. Dunque, anche se il metodo è convergente e il passo,
per esempio, k = 0.1 garantisce un errore locale proporzionale a k 2 = 0.01, il
metodo di Eulero non può essere usato con tale passo. Usando il metodo di
Eulero implicito sarebbe possibile invece usare un passo piccolo all’inizio e
poi, quando ormai la prima componente è trascurabile, si potrebbe incremen-
tare il passo, senza pericolo di esplosione della soluzione. Per questo semplice
problema, sarebbe possibile calcolare le due componenti separatamente. Nel
caso generale, però, il sistema non è disaccoppiato. Per l’analisi, ci si può
ricondurre, eventualmente in maniera approssimata, ad uno disaccoppiato e
ragionare per componenti. Infatti, se A è una matrice diagonalizzabile,
y ′ (t) = Ay(t) ⇔ z ′ (t) = Dz(t) ⇔ z(t) = exp(tD)z 0
ove AV = V D, D = diag{λ1 , λ2 , . . . , λd }, e y(t) = V z(t). Poi
• ℜ(λ1 ) ≪ ℜ(λ2 )
y ′ (t) = Ay(t)
Integratori esponenziali
I problemi di assoluta stabilità per semplici problemi lineari visti nel capitolo
precedente, portano alla ricerca di nuovi metodi. Consideriamo il sistema
differenziale (
y ′ (t) = Ay(t) + b, t > 0
y(t0 ) = y 0
La soluzione analitica è
y(t) = exp((t−t0 )A)y 0 +(t−t0 )ϕ1 ((t−t0 )A)b = y 0 +(t−t0 )ϕ1 ((t−t0 )A)(Ay 0 +b)
Infatti y(t0 ) = y 0 e
130
131
Dimostrazione. Si ha
Z tn+1
y n+1 = exp(kA)y n + exp((tn+1 − τ )A)g(tn )dτ
tn
ove si è posto g(t) = b(t, y(t)). Per la formula di variazioni delle costanti (7.2)
Z tn+1
y(tn+1 ) − exp(kA)y(tn ) − exp((tn+1 − τ )A)g(tn )dτ =
tn
Z tn+1
= exp(kA)y(tn ) + exp((tn+1 − τ )A)g(τ )dτ +
tn
Z tn+1
− exp(kA)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
= k ϕ2 (kA)g ′ (τn ) = O(k 2 )
2
ove
∂f
Jn = (y ), bn (y(t)) = f (y(t)) − Jn y(t)
∂y n
e applicarvi il metodo di Eulero esponenziale. Si arriva cosı̀ al metodo di
Eulero–Rosenbrock esponenziale
Esercizi
133
134 CAPITOLO 21. ESERCIZI
ϑ0
ϑ(t)
0
1 1
2 2
1 1
2
0 2
1 0 0 1
1 1 1 1
6 3 3 6
c A
bT
0
1 1
2 2
1 −1 2
0 1
1 2 1
6 3 6
PDEs
(Equazioni alle derivate
parziali)
136
Capitolo 22
Equazioni di
trasporto-diffusione-reazione
137
138 CAPITOLO 22. EQUAZIONI ADR
da cui
ψ ′ (t)
= −K (costante) ⇒ ψ(t) = Ae−Kt
ψ(t)
Per quanto riguarda φ(x), la soluzione generale è
√ √
−Kx
φ(x) = Be + Ce− −Kx
0 = φ(0) = B + C
√ √ ³ √ √ ´
−KL
0 = φ(L) = Be + Ce− −KL
= B e −KL − e− −KL
√ √
Se K ≤ 0, allora e −KL − e− −KL > 0 e dunque B = 0 (e anche C). Quindi
φ(x) = 0, ma in tal caso ψ(0)φ(x) 6= u0 (x). Se invece K = λ2 > 0, λ > 0,
allora ¡ ¢
φ(x) = B eiλx − e−iλx = 2Bi sin(λx) = B sin(λx)
(avendo ridefinito B) e poiché φ(L) = 0, l’unica possibilità non banale è
λ = jπ/L, j numero naturale non nullo. Pertanto, la funzione
µ 2 2 ¶ µ ¶
j π jπ
uj (t, x) = exp − 2 t sin x
L L
converge in [−L, L]. Poiché ū0 (x) è dispari, con riferimento al paragra-
fo 13.2.1,
Z L Z L µ ¶
−i 2πj(x + L)
u0m/2+1+j = ū0 (x)φm/2+1+j (x)dx = √ ū0 (x) sin dx =
−L 2L −L 2L
√ Z µ ¶
−i 2 L jπx
= √ ū0 (x) sin + jπ dx =
L 0 L
√ Z µ ¶
−i 2 L jπx
= √ u0 (x) sin + jπ dx
L 0 L
e
Z L Z L µ ¶
−i −2πj(x + L)
u0 m/2+1−j = ū0 (x)φm/2+1−j (x)dx = √ ū0 (x) sin dx =
−L 2L −L 2L
√ Z
i 2 L
= √ ū0 (x) sin(jπx/L + jπ)dx = −u0 m/2+1+j
L 0
da cui
+∞
X
ū0 (x) = u0 m/2+1+j φm/2+1+j (x) =
j=−∞
+∞
X cos(jπx/L + jπ) + i sin(jπx/L + jπ)
= u0m/2+1+j √ =
j=−∞
2L
−1
X cos(jπx/L + jπ) + i sin(jπx/L + jπ)
= u0m/2+1+j √ +
j=−∞
2L
+∞
X cos(jπx/L + jπ) + i sin(jπx/L + jπ)
+ u0 m/2+1+j √ =
j=1
2L
+∞
X cos(jπx/L + jπ) − i sin(jπx/L + jπ)
= −u0m/2+1+j √ +
j=1
2L
+∞
X cos(jπx/L + jπ) + i sin(jπx/L + jπ)
+ u0 m/2+1+j √ =
j=1
2L
∞ √
X 2
= u0 m/2+1+j √ i sin(jπx/L + jπ) =
j=1
L
∞ · Z µ ¶ ¸ µ ¶
X 2 L jπ jπ
= u0 (x) sin x dx sin x
j=1
L 0 L L
140 CAPITOLO 22. EQUAZIONI ADR
Si ha
Z L · ¸ Z L Z L 2
dE ∂ 1 2 ∂u ∂ u
= u (t, x) dx = u dx = u 2 dx
dt 0 ∂t 2 0 ∂t 0 ∂x
la soluzione approssimata, si ha
m
X m
X
û′j (t)φj (x) = ûj (t)λ2j φj (x)
j=1 j=1
ove λk = i(k − 1 − m/2)2π/(2L) (si deve lavorare infatti del dominio [−L, L],
dove si è prolungata per antisimmetria la funzione u0 (x)) √ e uˆ√
0 k sono i coeffi-
cienti di Fourier discreti di u0 prolungata (per inciso, si ha i 2 Luˆ0 m/2+1+k =
(−1)k ck ). Si trova, dunque,
2 π 2 t/L2
ûk (t) = e−(k−1−m/2) uˆ0 k , 1≤k≤m
142 CAPITOLO 22. EQUAZIONI ADR
da cui poi si ricostruisce û(t, x). Avevamo visto che la decomposizione di Fou-
rier si usa in caso di condizioni al bordo periodiche, mentre per l’equazione
del calore sono di Dirichlet nulle. Poiché però il dato iniziale è la funzione di-
spari ū0 (x), allora la soluzione ū(t, x) dell’equazione del calore nell’intervallo
[−L, L] è pure dispari. Infatti, posto v̄(t, x) = −ū(t, −x), si ha
∂v̄ ∂ ū
(t, x) = − (t, −x)
∂t ∂t
2
∂ v̄ ∂ 2 ū
(t, x) = (− · −) · − (t, −x)
∂x2 ∂x2
inoltre, v̄(t, −L) = v̄(t, L) = 0 e v̄(0, x) = −ū(0, −x) = ū0 (x). Dunque, pure
v̄(t, x) soddisfa l’equazione del calore. Ma questa è unica, quindi v̄(t, x) =
−ū(t, −x) = ū(t, x), cioè ū(t, x) è dispari. Quindi, ∂x ū(t, −x) = ∂x ū(t, x) e,
in particolare, ∂x ū(t, −L) = ∂x ū(t, L). Per quanto visto, la serie di Fourier di
ū(t, x) converge (i coefficienti uj (t) decadono a zero almeno come j 2 ) e, poiché
ogni troncata della serie è dispari e periodica, essa vale zero ai bordi x = −L
ed x = L (e, di conseguenza, anche in x = 0: questo fatto non è vero per
l’equazione originaria nel dominio [0, L], poiché lı̀ la soluzione non è dispari).
Dunque, si può usare il metodo di Fourier. Se però ū0 (x) non è periodica
(nel senso che non lo sono le derivate di ordine superiore al primo), allora
tale sarà la soluzione analitica e il metodo di Fourier non sarà spettralmente
convergente.
ove il termine g(u(t, x)) si chiama reazione e il termine s(t, x) sorgente, pre-
vede 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. Assumeremo
sempre che la condizione iniziale soddisfi le condizioni ai bordi. Vediamo
qualche esempio.
22.3. METODO DELLE LINEE 143
che vengono poi amplificati dal coefficiente 1/h2 . Dunque, con riferimento
alla condizione (19.2) per il metodo di Eulero, volendo usare questo metodo
per l’integrazione temporale occorrerebbe un passo temporale k minore di
(circa) h2 /2. Siccome il metodo di Eulero è del primo ordine, volendo che
l’integrazione temporale non sia meno accurata dell’approssimazione spaziale,
è giusto che il passo temporale sia proporzionale a h2 (cosı̀ che l’errore globale
sia O(k + h2 ) = O(h2 )). Per ridurre il numero di time steps, si può usare un
metodo di ordine più alto, per esempio un metodo di Runge–Kutta esplicito
di ordine 2. La restrizione sul time step è però la stessa (vedi la regione di
assoluta stabilità del metodo in Figura 19.3) del metodo di Eulero. Dunque,
ancora k dovrebbe essere proporzionale a h2 /2 (quindi il numero di time steps
non diminuisce) e l’errore globale è ancora O(k 2 + h2 ) = O(h4 + h2 ) = O(h2 ).
144 CAPITOLO 22. EQUAZIONI ADR
peclet.m s hema SG
soluzione di riferimento
FD entrate
s hema SG
soluzione di riferimento
FD entrate
1 1
0.5 0.5
u
u
0 0
-0.5 -0.5
0.9 0.91 0.92 0.93 0.94 0.95 0.96 0.97 0.98 0.99 1 0.9 0.91 0.92 0.93 0.94 0.95 0.96 0.97 0.98 0.99 1
x x
10−4
10−5
errore
10−6
10−7
10−8
101 102 103
passi temporali
22.5 Esercizi
1. Si calcoli la soluzione analitica dell’equazione del calore con sorgente
∂u ∂2u
(t, x) = 2
(t, x) + 2et sin(x), t > 0, x ∈ (0, π/2)
∂t ∂x
u(t, 0) = 0, t>0
∂u
(t, π/2) = 0, t>0
∂x
u(0, x) = sin(x), x ∈ (0, π/2)
usando differenze finite del secondo ordine nello spazio e il metodo dei
trapezi nel tempo. Si mostrino gli ordini spaziali e temporali della
convergenza alla soluzione analitica al tempo t∗ = 1.
• il metodo di Eulero
• il metodo di Eulero implicito
• il metodo dei trapezi
• il metodo di Heun
• il metodo Runge–Kutta di tableau in Tabella 21.1
reazione
∂u ∂u ∂ 2u
(t, x) + c (t, x) = d (t, x) + ρu(t, x)(u(t, x) − 1/2)(1 − u(t, x)), t > 0, x ∈ (0, 1)
∂t ∂x ∂x2
u(t, 0) = 0, t>0
∂u
(t, 1) = 0, t>0
∂x
u(0, x) = 5x(1 − x)2 , x ∈ (0, 1)
Temi d’esame
154
155
29-06-2010
0
1 1
2 2
1 1
2
0 2
1 0 0 1
3 5 7 13 1
4 32 32 32
− 32
− 12 7
3
7
3
13
6
− 16
3
u(0, x) = 4x(1 − x), x ∈ (0, 1)
u(t, 0) = u(t, 1) = 0, t>0
20-07-2010
156
15-09-2010
di tableau √ √
3+ 3 3+ 3
c1 a1,1 a1,2 6√ 6√ √
0
3− 3
c2 a2,1 a2,2 = 6
− 33 3+ 3
6
b1 b2 1 1
2 2
per la soluzione di un problema ai valori iniziali
(
y ′ (t) = f (t, y(t)), t ≥ t0
y(t0 ) = y 0
24-02-2011
14. Si risolva mediante un metodo di shooting il problema ai limiti
′′ 1
u (x) = (32 + 2x3 − u(x)u′ (x)), x ∈ (1, 3)
8
u(1) = 17
43
u(3) =
3
Sapendo che la soluzione analitica è u(x) = x2 + 16/x, si determini
sperimentalmente il numero minimo di passi temporali per avere un
errore in norma infinito minore di 10−2 .
15. Si applichi il metodo delle linee al problema di diffusione non omogeneo
∂u ∂2u 5 x
(t, x) = 2
(t, x) + et sin , t > 0, x ∈ (0, π)
∂t ∂x 4 2
x
u(0, x) = sin , x ∈ (0, π)
2
u(t, 0) = 0, t≥0
∂u (t, π) = 0,
t≥0
∂x
e si mostri l’ordine di convergenza temporale del metodo esponenziale
punto medio con cui viene approssimata la soluzione analitica u(t∗ , x)
al tempo t∗ = 1.
21-06-2011
16. Si risolva il seguente problema differenziale
′′
u (x) + u(x) = 2 cos(x), x ∈ (0, 1]
u(0) = 0
u′ (0) = 0
11-07-2011
15-09-2011
y(0) = 1
′
y (0) = 0
24-02-2012
0 0
1 1
1 2 2
1 1
2 2
21-06-2012
18 9 2 6
y n+3 − y + y − y n = k f (tn+3 , y n+3 )
11 n+2 11 n+1 11 11
18 9 3 6
y n+3 − y n+2 + y − y n = k f (tn+3 , y n+3 )
11 11 n+1 11 11
si dica quale è consistente e perché e se ne determini numericamente
l’ordine.
05-07-2012
164
10-09-2012
usando differenze finite del secondo ordine con passo h = 1/50 nello spa-
zio e il metodo di Eulero esplicito nel tempo fino a t∗ = 1. Perché serve
un passo temporale minore di circa 1/20000 per avere convergenza? Si
mostri infine il corretto ordine di convergenza temporale.
24-09-2012
Bibliografia
166
Bibliografia
[4] E. Hairer and G. Wanner, Solving Ordinary Differential Equations II, Stiff
and Differential-Algebraic Problems, Springer, Second Revised Edition,
2002.
167