Analisi Numerica (A.A. 2014-2015) - Università di Roma · Analisi Numerica (A.A. 2014-2015)...
112
Analisi Numerica (A.A. 2014-2015) Appunti delle lezioni: Equazioni differenziali ordinarie
Transcript of Analisi Numerica (A.A. 2014-2015) - Università di Roma · Analisi Numerica (A.A. 2014-2015)...
Analisi Numerica
(A.A. 2014-2015)
Appunti delle lezioni:
Equazioni differenziali ordinarie
Docente Vittoria Bruni
Email: [email protected]
Ufficio: Via A. Scarpa,
Pal. B, I piano, Stanza n. 16
Tel. 06 49766648
Ricevimento: Martedı 9.30-10.30
Testi consigliati:
Calcolo Numerico, L. Gori, Ed. Kappa, 2006
Esercizi di Calcolo Numerico, L. Gori-M.L. Lo Cascio, F. Pitolli, Ed. Kappa, 2007
Il materiale didattico e disponibile sul sito
http://www.sbai.uniroma1.it/users/bruni-vittoria
nella pagina dedicata al corso Analisi Numerica
1
Equazioni differenziali
Le equazioni differenziali rappresentano uno strumento efficace per lacostruzione di modelli matematici per la simulazione di sistemi dinamiciin diversi ambiti dell’ingegneria ma anche della biologia, della chimica,della fisica, dell’informatica e delle scienze sociali.
Sono equazioni che esprimono un legame tra una o piu funzioni incog-nite e le loro derivate.
Se le derivate si riferiscono ad una sola variable, allora parliamo diequazioni differenziali ordinarie (spesso indicate con l’acronimo ODE),per esempio
y′(t) = ay(t),
con y funzione incognita, t variabile indipendente e a costante; se, in-vece, le derivate si riferiscono a piu variabili, allora parliamo di equazionidifferenziali alle derivare parziali (in genere indicate con l’acronimoPDE), per esempio
ut − auxx = 0,
con u funzione incognita, t e x variabili indipendenti e a costante.
2
Eq. differenziali ordinarie: modello matematico
Il moto di una particella di massa m attaccata all’estremita di una
molla di costante elastica k e descritto dall’equazione differenziale
lineare del secondo ordine (equazione dell’oscillatore armonico
semplice smorzato):
md2x
dt2+ b
dx
dt+ kx = 0
−b dxdt
: forza di attrito
−k x: legge di Hooke
Dall’analisi matematica sappiamo che la soluzione ”esatta” e
x(t) = xm e−bt/2m cos(ωm t+ ϕm) ωm =
√k
m− b2
4m2:
pulsazione
dell’oscillatore
L’ampiezza xm e la fase ϕm dell’oscillazione sono individuate dalle
condizioni iniziali.
3
Eq. differenziali ordinarie: problema di Cauchy
• Un problema reale richiede l’assegnazione delle condizioni ini-
ziali.
Problema di Cauchy:
md2x
dt2+ b
dx
dt+ kx = 0
x(0) = x0 c.i.
dx
dt(0) = v0 c.i.
⇒
x(t) = xm e−bt/2m cos(ωm t+ ϕm)
xm =x0
cosϕm
tanϕm =v0
x0
2m
b
Esempio:
x(0) = 1
dx
dt(0) = 0
0 5 10 15 20−2
−1.5
−1
−0.5
0
0.5
1
1.5
t
dx/dtx
4
Eq. differenziali ordinarie: modello matematicoUn modello analogo descrive un problema di origine diversa: determinare la caricaQ(t) di un consendatore C nel circuito chiuso illustrato in figura, in cui R e laresistenza, L l’induttanza, E una sorgente di forza elettromotrice e S l’interruttore
In questo caso si applica la seconda legge di Kirchoff che stabilisce che in un circuitochiuso la differenza di potenziale fornita attraverso la sorgente E e uguale alla sommadelle differenze di potenziale misurate nel resto del circuito, ovvero
E(t) = RI(t) + LdI
dt(t) +
Q(t)
C
con I l’intensita di corrente. Poiche I(t) = dQdt
(t), si ha
Ld2Q
dt2+R
dQ
dt+Q
C= E(t)
5
Eq. differenziali ordinarie: problema di Cauchy
• Nei problemi reali l’espressione dell’equazione differenziale e
complicata e, in genere, non si riesce a calcolare esplicitamente
la soluzione.
Esempio
Le oscillazioni di un pendolo possono essere descritte dall’equazione
differenziale del secondo ordine non lineare
d2θ
dt2− g
Lsin θ = 0
dove L e la lunghezza del pendolo, g e l’accelerazione di gravita e θ e
l’angolo tra il pendolo e la verticale. Il problema e completato con le
condizioni iniziali
θ(t0) = θ0 θ′(t0) = v0
Dati: g = 9.81 m/s2, L = 0.6 m, θ0 = π/6, v0 = 0
6
Sistemi di equazioni differenziali: esempio
Il modello preda-predatore di Lotka-Volterra e il piu semplice modello
che descrive la competizione tra due specie. Consiste in una coppia di
equazioni differenziali in cui y1(t) rappresenta il numero di prede e y2(t)
rappresenta il numero di predatori. L’ interazione tra le due specie e
descritta da un termine proporzionale a entrambe le popolazioni:y′1(t) = k1
(1− y2(t)
µ2
)y1(t), t > 0
y′2(t) = −k2
(1− y1(t)
µ1
)y2(t)
dove k1, k2, µ1, µ2 sono costanti positive e y1,0, y2,0 rappresentano
le condizioni iniziali.
In assenza di predatori le prede crescono in modo esponenziale; i preda-
tori invece, in assenza di prede, muoiono rapidamente.
Oss. Se y1,0 = µ1, y2,0 = µ2, allora y′1 = 0 e y′2 = 0, ovvero (µ1, µ2)
e un punto di equilibrio. Per particolari valori di µ1, µ2 la soluzione e
periodica.
7
Equazioni differenziali di ordine n
y(n)(x) = f(x, y(x), y′(x), y′′(x), . . . , y(n−1)(x))
y(k)(x0) = yk k = 0,1, . . . , n− 1 Condizioni iniziali
Un’equazione differenziale di ordine n puo essere ricondotta a un
sistema di n equazioni differenziali del primo ordine.
y(n)(x) = f(x, y(x)︸ ︷︷ ︸y1(x)
, y′(x)︸ ︷︷ ︸y2(x)
, y′′(x)︸ ︷︷ ︸y3(x)
, . . . , y(n−1)(x)︸ ︷︷ ︸yn(x)
)
y1(x) = y(x)
y′1(x) = y2(x)y′2(x) = y3(x).........................y′n−1(x) = yn(x)y′n(x) = f(x, y1(x), y2(x), y3(x), . . . , yn(x))
y10 = y0y20 = y1..................yn0 = yn−1
Condizioni
iniziali
8
Eq. differenziali ordinarie del primo ordine /1
Problema di Cauchy:
{y′(x) = f(x, y(x))y(x0) = y0 condizione iniziale
Esempio: Modello per descrivere la crescita di una popolazione
f(x, y(x)) = Ky(x) K > 0 −→{y′(x) = Ky(x)y(x0) = y0
Soluzione esatta:
y(x) = y0eK(x−x0)
xx0
y(x)
y0
9
Eq. differenziali ordinarie del primo ordine /2
Nei problemi applicativi, l’espressione di f(x, y(x)) e complicata equindi non si puo calcolare esplicitamente la soluzione del problemadi Cauchy oppure la sua espressione e data tramite funzioni non ele-mentari.
−→ si deve ricorrere a un metodo numerico per approssimare y(x).
Esempio
Sia P (t) il numero di individui in una popolazione al tempo t, misuratoin anni. Se il tasso medio di nascita b e costante e il tasso di morte d eproporzionale al numero di individui presenti al tempo t, allora il tassodi crescita e dato dall’equazione differenziale del primo ordine
dP (t)
dt= bP (t)− d(t)P (t), d(t) = kP (t),
chiamata equazione logistica.
Assumendo P (0) = 50976, b = 2.9× 10−2 e k = 1.4× 10−7, calcolareil numero di individui dopo 5 anni.
10
Esistenza e unicita dellasoluzione del problema di Cauchy
Prima di risolvere numericamente un’equazione differenziale bisognaessere sicuri che il problema di Cauchy ammetta un’unica soluzione
y(x).
Definizione. Una funzione f(x, y) si dice lipschitziana in y uniforme-
mente rispetto a x in D ⊂ IR2, se esiste una costante
L > 0 tale che
|f(x, y1)− f(x, y2)| ≤ L |y1 − y2| ∀(x, y1), (x, y2) ∈ D.
Nota. Una condizione sufficiente e che fy esista e sia limitata in D.
Esempio: f(x, y) = K y → fy(x, y) = K → e limitata
→ f e lipschitziana
Al contrario, la funzione f(x) = |x| e continua ma non C1 in un intervallo I limitato e
simmetrico rispetto allo 0, in quanto la derivata non e continua in x = 0. Tuttavia, la
funzione e Lipschitziana con L = 1 in quanto ∀ x1, x2 ∈ I, si ha ||x1|−|x2|| ≤ |x1−x2|.11
Teorema 1. Sia f definita e continua in Q:
Q :={
(x, y) ∈ IR2 |x ∈ [x0, x0 + a]; y ∈ [y0 − b, y0 + b]}
e sia lipschitziana in y, uniformemente
rispetto a x
x0
y(x)
x0+a
y0−b
y0+b
y0
Q
⇒ esiste un’unica soluzione y(x) ∈ C1[x0, x0 + β] del problema di
Cauchy in I = [x0, x0 + β] con β = min[a,
b
M
], M = maxQ|f(x, y)|.
Teorema 2. Sia f definita e continua in S:
S :={
(x, y) ∈ IR2 |x ∈ I = [x0, x0 + β]; y ∈ IR}
e sia lipschitziana in y, uniformemente
rispetto a x
x0
y(x)
x0+β
y0
S
⇒ esiste un’unica soluzione y(x) ∈ C1[x0, x0 + β] del problema di
Cauchy in I = [x0, x0 + β].
12
Ben-posizione del problema di Cauchy
Per poter risolvere numericamente il problema di Cauchy e necessarioanche che il problema sia ben posto.
Definizione. Il problema di Cauchy
{y′(x) = f(x, y(x))y(x0) = y0 condizione iniziale
e ben posto se, dette y(x; y0) e y(x; y0 +δ) le soluzioni con condizioni
iniziali y(x0) = y0 e y(x0) = y0 + δ rispettivamente, si ha
|y(x; y0)− y(x; y0 + δ)| < ε x ∈ [a, b]
dove ε > 0 e una prefissata tolleranza, purche δ = δ(ε) sia sufficiente-mente piccolo.
Nota. Se f soddisfa le condizioni del Teorema 2, il problema di
Cauchy e ben posto.
13
Esempi
• Il seguente problema di Cauchy:{y′(t) = −y + t+ 1 t ∈ [0,10]y(0) = 1
ammette soluzione unica in I = [0,10]. Infatti, f(t, y) = −y + t+ 1
e definita e continua in [0,10]×R Inoltre, fy(t, y) = −1, ∀ t ⇒|fy| = 1, ∀ t ∈ I ∀ y ∈ R
• Il seguente problema di Cauchy:{y′(t) = 1 + tsin(ty) t ∈ [0,2]y(0) = 0
ammette soluzione unica in I = [0,2]. Infatti, f(t, y) = 1+tsin(ty) e
definita e continua in [0,2]×R Inoltre, fy(t, y) = t2cos(ty), ∀ t ⇒|fy| = |t2cos(ty)| ≤ t2 ≤ 22 = 4, ∀ t ∈ I ∀ y ∈ R
14
Esempio
Si consideri il seguente problema di Cauchy{y′(t) = 3
2y13
y(0) = 0
Le soluzioni del problema sono:
y(t) = 0 e y(t) = t32
Inoltre, per a ∈ R, a > 0, la funzione
ya(t) =
{(t− a)
32 t > a
0 t = 0
e soluzione del problema di Cauchy.
La soluzione, dunque, non e unica. Infatti, f(t, y) non e lipshitziana
rispetto a y.
15
Metodo di Eulero
Problema di Cauchy:
{y′(x) = f(x, y(x)) x ∈ [x0, x0 + β]y(x0) = y0 condizione iniziale
Discretizzazione di I: xi = x0 + ih i = 0, . . . , n h =β
n
Metodo numerico: I valori esatti y(xi) vengono approssimati con i
valori yi.
Sviluppo in serie di Taylor:
y(x1) = y(x0 + h) = y(x0) + y′(x0)h+1
2y′′(τ1)h2 =
= y(x0) + f(x0, y(x0))h+1
2y′′(τ1)h2
τ1 ∈ [x0, x1]
Soluzione approssimata:
y1 = y0 + hf(x0, y0)
Errore globale di troncamento:
e1(x1) = y(x1)− y1 = P1T1
x
y
P0
P1
P2
P3
T1
T2
T3
x0 x
0+h x
0+2h x
0+3h
y0(x
0)
y1
y2
y3
y(x0+h)
y(x0+2h)
y(x0+3h)
t0
r1
r2
16
Algoritmo del metodo di Eulero
Algoritmo:
yi+1 = yi + hf(xi, yi) i = 0,1, . . . , n
Errore globale di troncamento:
ei = y(xi)− yi = PiTi
x
y
P0
P1
P2
P3
T1
T2
T3
x0 x
0+h x
0+2h x
0+3h
y0(x
0)
y1
y2
y3
y(x0+h)
y(x0+2h)
y(x0+3h)
t0
r1
r2
L’errore globale di troncamento ha due contributi:
• l’errore locale di troncamento, dovuto al fatto che la soluzione
esatta y(x) viene approssimata localmente con una retta:
R(xi, y(xi);h; f) =1
2h2y′′(τi) τi ∈ [xi, xi+1]
R(xi, y(xi);h; f) = O(h2) → Primo ordine
• l’accumularsi degli errori locali di troncamento, per cui al generico
passo i ≥ 1 ci si muove lungo la retta ri che e una approssimazione
della retta tangente alla soluzione in Pi = [xi, y(xi)].17
Convergenza del metodo di Eulero
Definizione. Un metodo numerico per la soluzione approssimata diun’equazione differenziale e convergente se
limh→0
max1≤i≤n|ei| = 0 ⇔ lim
i→∞ yi = y(x) con x = x0 + ih
Dalla figura e evidente che, se si riduce
il passo h, si riduce anche l’errore
ei = PiTi. Per di piu
limh→0
max1≤i≤n|ei| = 0
⇒ Il metodo di Eulero e convergentex
y
P0
P2
P4
T1
T2
T3
x0
x0+h x
0+2h
y(x0)
y1
y3
y4
y(x0+h/2)
y(x)
t0
r1
P1
r2
r3 T
4
P3
x0+h/2 x
0+3h/2
y2
y(x0+h)
r4
Un metodo numerico e convergente se e
• Consistente:
limh→0
R(x, y;h; f)
h= 0
+
• Stabile:
l’accumularsi degli errori locali di
troncamento si mantiene limitato18
Esercizio
Dato il problema di Cauchy (PC){
y′ = |y − x|, x ∈ [0,2]y(0) = 2
verificare l’ esistenza e unicita della soluzione in S = [0,2] × R e ap-
prossimarla con il metodo di Eulero con passi h1 = 0.1 e h2 = 0.05 per
x ∈ [0,1]
19
Soluzione
Per verificare se vale il teorema di esistenza e unicita della soluzione
bisogna verificare se f(x, y) = |y − x| e lipschitziana in S = [0,2] × R,
cioe se esiste una costante L > 0 (costante di Lipschitz) tale che
|f(x, y1)− f(x, y2)| ≤ L|y1 − y2|, ∀ (x, y1), (x, y2) ∈ S.
20
Per f(x, y) = |y − x| si ha
f(x, y) =
{y − x, x < yx− y, x ≥ y
⇒ fy(x, y) =
{1, x < y−1, x > y
∀ (x, y) ∈ S = [0,2]×R.
Quindi, |fy| = 1 x 6= y e |y − x| = 0 x = y (la derivata di f non e
definita in y = x). Poiche fy risulta limitata in x 6= y, si puo concludere
che f(x, y) e lipschitziana in S con costante di Lipschitz pari a 1. Dalla
lipschitzianita di f segue anche che il metodo di Eulero e convergente
per x ∈ [0,2].
21
Scegliendo il passo h = 0.1, definiamo i punti equidistanti
xi = x0 + ih
in cui si approssimeranno i valori della soluzione del PC, cioe
yi = yi−1 + hf(xi−1, yi−1).
Poiche si chiede la soluzione nell’intervallo [0,1], sara necessario ese-
guire
n =b− ah
=1
0.1= 10
passi del procedimento suddetto.
22
y1 = y0 + hf(x0, y0) = 2 + 0.1|2− 0| = 2.2, x1 = x0 + h = 0.1
y2 = y1 + hf(x1, y1) = 2.2 + 0.1|2.2− 0.1| = 2.41, x2 = x0 + 2h = 0.2
y3 = y2 + hf(x2, y2) = 2.41 + 0.1|2.41− 0.2| = 2.631, x3 = x0 + 3h = 0.3
y4 = y3 + hf(x3, y3) = 2.631 + 0.1|2.631− 0.3| = 2.8641, x4 = x0 + 4h = 0.4
y5 = y4 + hf(x4, y4) = 2.8641 + 0.1|2.8641− 0.4| = 3.11051, x5 = x0 + 5h = 0.5
...
...
y10 = y9 + hf(x9, y9) = 4.257947691 + 0.1|4.257947691− 0.9| = 4.5937424601,
x10 = x0 + 10h = 1
23
Per valutare l’accuratezza della approssimazione prodotta, confrontia-
mo i valori ottenuti con la soluzione esatta y∗(x) = ex+x+1 valutando
l’errore globale ei = |yi − y∗i |, cioe
xi yi y∗i ei0.0 2.0000000000 2.00000000000000 0.000000000000000.1 2.2000000000 2.20517091807565 0.005170918075650.2 2.4100000000 2.42140275816017 0.011402758160170.3 2.6310000000 2.64985880757600 0.018858807576000.4 2.8641000000 2.89182469764127 0.027724697641270.5 3.1105100000 3.14872127070013 0.038211270700130.6 3.3715610000 3.42211880039051 0.050557800390510.7 3.6487171000 3.71375270747048 0.065035607470480.8 3.9435888100 4.02554092849247 0.081952118492470.9 4.2579476910 4.35960311115695 0.101655420156951.0 4.5937424601 4.71828182845905 0.12453936835904
24
L’errore di stima cresce man mano che xi si avvicina a 1.
25
Consideriamo ora il passo di discretizzazione h = 0.05. In questo caso
eseguiamo 20 passi del metodo di Eulero:
xi yi y∗i ei0.00 2.00000000000000 2.00000000000000 00.05 2.10000000000000 2.10127109637602 0.001271096376020.10 2.20250000000000 2.20517091807565 0.002670918075650.15 2.30762500000000 2.31183424272828 0.004209242728280.20 2.41550625000000 2.42140275816017 0.005896508160170.25 2.52628156250000 2.53402541668774 0.007743854187740.30 2.64009564062500 2.64985880757600 0.009763166951000.35 2.75710042265625 2.76906754859326 0.011967125937010.40 2.87745544378906 2.89182469764127 0.014369253852210.45 3.00132821597852 3.01831218549017 0.016983969511650.50 3.12889462677744 3.14872127070013 0.019826643922690.55 3.26033935811631 3.28325301786740 0.022913659751080.60 3.39585632602213 3.42211880039051 0.026262474368380.65 3.53564914232324 3.56554082901390 0.029891686690660.70 3.67993159943940 3.71375270747048 0.033821108031080.75 3.82892817941137 3.86700001661267 0.038071837201310.80 3.98287458838194 4.02554092849247 0.042666340110530.85 4.14201831780103 4.18964685192599 0.047628534124960.90 4.30661923369108 4.35960311115695 0.052983877465870.95 4.47695019537564 4.53570965931585 0.058759463940211.00 4.65329770514442 4.71828182845905 0.06498412331463
26
Riducendo il passo di discretizzazione, l’errore si e ridotto di circa
la meta
Infatti, l’errore locale di troncamento e
R(xi, y(xi);h; f) =1
2h2y
′′(ηi) ηi ∈ [xi, xi+1]
cioe
R(xi, y(xi);h; f) =1
2h2eηi
27
Esercizio
Scrivere la funzione matlab eulero.m che implementi il metodo di
Eulero esplicito. La funzione ha come parametri di input il punto
iniziale (x0, y0) la funzione f(x, y(x)), il passo di discretizzazione h e il
numero di iterazioni da eseguire n step. L’output e una matrice T le cui
righe contengono rispettivamente i nodi e l’approssimazione prodotta
dal metodo nei nodi.
Usare la funzione per risolvere l’esercizio precedente usando i passi
h = 0.1, 0.05, 0.025, 0.01. Si calcolino gli errori commessi con-
frontando l’approssimazione prodotta con la soluzione esatta e verifi-
care che l’errore massimo nell’intervallo di stima e una funzione lineare
rispetto ad h.
28
function [T]=eulero(x0,y0,fun,h,n_step)
%
% [T]= eulero(x0,y0,fun,h,n_step)
% Soluzione di un’ equazione differenziale
% del primo ordine con il metodo di Eulero
%
% Dati di input:
% x0, y0: condizione iniziale
% n: numero dei passi da eseguire
% h: passo di discretizzazione
% fun: espressione di f(x,y)
%
% OUTPUT
% T = matrice di dimensione 2xn_step
% la prima riga contiene il vettore dei nodi,
% la seconda riga contiene il vettore delle approssimazioni
% della soluzione nei nodi
%
29
xi(1) = x0;
yi(1) = y0;
for i = 2:n_step+1
yi(i) = yi(i-1) + h*fun(xi(i-1),yi(i-1));
xi(i) = x0 + (i-1)*h;
end
T = [xi;yi];
30
Script matlab
y_true = @(x)[exp(x)+x+1];x0 = 0;y0 = 2;fun = @(x,y)[abs(y-x)];intervallo = 1;h = [0.1 0.05 0.025 0.01];color = [’b’,’r’,’g’,’k’];EH = [];figure,for i = 1:length(h)
n_step = 1/h(i);[T]=eulero(x0,y0,fun,h(i),n_step);err = abs(T(2,:)-y_true(T(1,:)));EH = [EH max(err)];hold on,plot(T(1,:),err,color(i))xlabel(’x_i’)ylabel(’ERRORE’)
endlegend(’0.1’,’0.05’,’0.025’,’0.01’)figure, loglog(h,EH)xlabel(’h’)title(’grafico di e_h’)
31
Grafico dell’errore di approssimazione del metodo di Eulero al variare
del passo di discretizzazione h.
32
Grafico del massimo dell’errore di approssimazione al variare del passo
di discretizzazione h e confronto delle soluzioni prodotte.
33
Metodi di ordine superiore al primo
Sono basati sull’uso dello sviluppo in serie di Taylor di y(x) arrestato
ad un certo ordine
Sviluppo in serie di Taylor di ordine m (y ∈ Cm+1[x0, xn])
y(xi+1) = y(xi) + h
m∑
k=1
y(k)(xi)
k!hk−1 +
hm+1
(m+ 1)!y(m+1)(τi) τi ∈ (xi, xi+1)
Schema di Taylor di ordine m
Trascurando l’errore e sostituendo i valori approssimati si ottiene l’algoritmo
yi+1 = yi + hm∑
k=1
f(k−1)(xi, yi)
k!hk−1
34
Errore locale di troncamento:
R(xi, y(xi);h; f) =hm+1
(m+ 1)!y(m+1)(τi)
Le derivate della f possono essere calcolate ricorsivamente da{f(k)(x, y(x)) = f
(k−1)x (x, y(x)) + f
(k−1)y (x, y(x)) f(x, y(x)) k ≥ 1
f(0)(x, y(x)) := f(x, y(x))
35
Metodi di ordine superiore al primo: Runge-Kutta
Si ottengono muovendosi lungo una retta che ha come coefficiente
angolare una combinazione lineare di valori assunti da f(x, y) in punti
opportuni dell’intervallo [xi, xi + h].
Metodi di Runge-Kutta del secondo ordine
Relazione esatta
y(x+h) = y(x)+hΦ(x, y;h; f)+R(x, y;h; f) → R(x, y;h; f) = O(h2+1)
con
Φ(x, y;h; f) = a1 k1(x, y) + a2 k2(x, y)
k1(x, y) = f(x, y)k2(x, y) = f (x+ λh, y + λh k1(x, y))
36
Sviluppo in serie di Taylor
y(x+ h) = y(x) + y′(x)h+1
2y′′(x)h2 + 0(h3) =
= y(x) + f(x, y)h+1
2
(fx(x, y) + fy(x, y) f(x, y)
)h2
︸ ︷︷ ︸hΦ(x;y;h;f)
+O(h3)
Si osserva che la funzione incremento cosı definita dipende anche da
fx e fy. Tuttavia, poiche f(x, y) + fx(x, y)λh + fy(x, y)λhf(x, y) e lo
sviluppo in serie di Taylor attorno al punto (x, y) della funzione
f(x+ λh, y + λhf(x, y)),
si ha
k2(x, y) = f(x, y) + fx(x, y)λh+ fy(x, y)λh f(x, y) +O(h2)
37
R(x, y;h; f) = y(x+ h)− y(x)− hΦ(x, y;h; f) =
= h f(x, y) +h2
2
(fx(x, y) + fy(x, y) f(x, y)
)+O(h3)−
−h(a1f(x, y)+
a2
(f(x, y) + fx(x, y)λh+ fy(x, y)λh f(x, y) +O(h2)
))=
= hf(x, y)[1− a1 − a2] + h2fx(x, y)(
12 − a2λ
)+
+h2fy(x, y)f(x, y)(
12 − a2λ
)+O(h3)
Affinche il metodo sia del secondo ordine devono essere nulli itermini in h e h2
1− a1 − a2 = 0
1
2− a2 λ = 0
che e un sistema di 2 equazioni in 2 incognite ⇒ ammette infinite soluzioni
38
Esempio: metodo di Heun
Il metodo di Heun si ottiene per una particolare scelta dei tre para-metri a1, a2, λ
a1 = a2 =1
2λ = 1
yi+1 = yi +h
2
(f(xi, yi) + f(xi + h, yi + hf(xi, yi))
)i = 0,1, . . . , n
y0 = y(x0)
Errore locale di trocamento:
R(x, y;h; f) = O(h3)
⇒ Secondo ordine
x
y
P1
T1
y(x)
t0
r0
t1
P0 y
0
x0
y(x1)
y1
39
Esempio: metodo di Eulero modificato
Il metodo di Eulero modificato si ottiene per un’altra particolare
scelta dei tre parametri a1, a2, λ
a1 = 0 a2 = 1 λ =1
2
yi+1 = yi + h
(f(xi +
h
2, yi +
h
2f(xi, yi))
)i = 0,1, . . . , n
y0 = y(x0)
40
Metodo di Runge-Kutta del 4o ordine
yi+1 = yi +
h
6(k1 + 2k2 + 2k3 + k4) i = 0,1, . . . , n
y0 = y(x0)
k1 = f(xi, yi)
k2 = f
(xi +
h
2, yi +
h
2k1
)
k3 = f
(xi +
h
2, yi +
h
2k2
)
k4 = f(xi + h, yi + hk3)
Errore locale di trocamento:
R(x, y;h; f) = O(h5) ⇒ Quarto ordine
41
Metodo di Runge-Kutta a r stadi
Piu in generale e possibile costruire metodi di Runge Kutta a r stadi
yi+1 = yi + h
r∑
l=1
alkl i = 0,1, . . . , n
y0 = y(x0)
con
k1 = f(xi, yi)
kl = f
xi + λih, yi + h
l−1∑
j=0
bijkk
che risultano consistenti ser∑
l=1
al = 1.
Inoltre, l’ordine del metodo e
p = r, r = 1, ...,4p = r − 1, r = 5, ...,7p = r − 2, r ≥ 8
42
Esempio
Problema di Cauchy:
{y′(x) = y(x)− x x ∈ [0,2]y(0) = 2
→ Soluzione esatta: y(x) = ex+x+1
0 0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 1.8 22
3
4
5
6
7
8
9
10
11
x
y
Metodo di Eulero
Metodo di Heun
Metodo di R−K classico
Grafico della soluzione
0 0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 1.8 20
0.2
0.4
0.6
0.8
1
1.2
1.4
Metodo di R−K classico
Metodo di Heun
Metodo di Eulero
x
Grafico dell’errore
43
Metodi one-step espliciti
I metodi di Eulero, Heun e Runge-Kutta classico sono tutti metodi
one-step espliciti, cioe metodi in cui per il calcolo di yi+1 si utilizza
solo il valore approssimato yi:
Metodi one-step espliciti:
yi+1 = yi + hΦ(xi, yi;h; f)︸ ︷︷ ︸ i = 0,1, . . . , n
Funzione incrementoy0 = y(x0)
Metodo di Eulero: Φ(x, y;h; f) = f(x, y(x))
Metodo di Heun: Φ(x, y;h; f) =1
2[f(x, y(x)) + f(x+ h, y(x) + hf(x, y(x))]
44
Errore locale di troncamento
L’errore locale di troncamento e definito come la differenza tra lasoluzione esatta y(x + h) e l’approssimazione fornita dallo schemanumerico quando questo utilizza i valori esatti y(x):
R(x, y(x);h; f) = y(x+ h)−(y(x) + hΦ(x, y(x);h; f)
)
Rappresenta l’errore dovuto al fatto di aver approssimato localmentela soluzione esatta con una retta.
Errore locale unitario di troncamento:
τ(x, y(x);h; f) =R(x, y(x);h; f)
h=y(x+ h)− y(x)
h−Φ(x, y(x);h; f)
Ordine del metodo: permette di valutare quanto e accurata l’approssi-mazione ottenuta. E’ definito come il massimo intero p tale che
R(x, y(x);h; f) = O(hp+1) ⇔ τ(x, y(x);h; f) = O(hp)
ovvero ∃ Cp costante : Cphp+1
• Un metodo di ordine p e esatto per tutti i polinomi fino al grado p.
45
Consistenza dei metodi one-step espliciti
Consistenza: limh→0
τ(x, y(x);h; f) = limh→0
R(x, y;h; f)
h= 0
I metodi one-step espliciti si ottengono discretizzando una relazione
esatta del tipo:
y(x+ h) = y(x) + hΦ(x, y(x);h; f) +R(x, y(x);h; f)︸ ︷︷ ︸↓
Errore locale di troncamento
limh→0
τ(x, y;h; f) = limh→0
[y(x+ h)− y(x)
h−Φ(x, y(x);h; f)
]
= y′(x)− limh→0
Φ(x, y(x);h; f) = 0
⇒ Un metodo one-step esplicito e consistente se e solo se
limh→0
Φ(x, y(x);h; f) = y′(x) = f(x, y(x))
46
In particolare,
• il metodo di Eulero e consistente. Infatti,
limh→0
Φ(x, y(x);h; f) = limh→0
f(y, y(x)) = f(x, y(x))
• il metodo di Heun e consistente.
Infatti,
limh→0
Φ(x, y(x);h; f) = limh→0
1
2(f(y, y(x))+f(x+h, y(x)+h)) = f(x, y(x))
47
Stabilita dei metodi one-step espliciti
Errore globale di troncamento: ei = y(xi)− yiConvergenza: lim
h→0max
0≤i≤n|ei| = 0 ⇔ limi→∞ yi = y(x) x = x0 + ih
ei+1 = y(xi+1)− yi+1 == (y(xi) + hΦ(xi, y(xi);h; f) +R(xi, y(xi);h; f))−− (yi + hΦ(xi, yi;h; f))
• Dalla definizione di ordine di un metodo si ha:
R(xi, y(xi);h; f) = O(hp+1)⇒ |R(xi, y(xi);h; f)| ≤ Cphp+1
Esempio: per il metodo di Eulero |R| = |12h2y′′(x)| ≤ M
2h2, dove M e il
massimo di |y′′(x)| nell’intervallo di integrazione ⇒ C1 = M/2
• Poiche Φ e una combinazione lineare di valori di f , dalla lip-shitzianita di f segue la lipshitzianita di Φ:
|Φ(x, y1)−Φ(x, y2)| ≤ L|y1 − y2|
48
Utilizzando le maggiorazioni trovate e trascurando gli errori di arro-tondamento si ha:
ei+1 = (y(xi)− yi)︸ ︷︷ ︸ei
+h (Φ(xi, y(xi);h; f)−Φ(xi, yi;h; f)︸ ︷︷ ︸≤L|y(xi)−yi|
+O(hp+1)
⇒{|ei+1| ≤ |ei|(1 + hL) + Cphp+1 i = 0,1, ...e0 = y(x0)− y0 = 0
Si associa alla disuguaglianza ottenuta una equazione alle differenzedel primo ordine:{ti+1 = ti(1 + hL) + Cphp+1 i = 0,1, ...t0 = 0
⇒ ti = Cphp(1 + hL)i − 1
L
Poiche |ei| ≤ ti, i = 0,1, ..., e (1 + α)i < eiα per α > 0, si deduce chel’errore globale di troncamento ha la limitazione
|ei| ≤ Cp hpeL(xi−x0) − 1
LL: costante di Lipschitz di f(x, y) Cp: costante dipendente dal metodo one-step.
49
Convergenza dei metodi one-step espliciti
Un generico metodo per risolvere numericamente un’equazione dif-
ferenziale e convergente se e consistente (l’errore locale unitario
di troncamento tende a zero quando h tende a zero) e stabile (la
propagazione degli errori durante lo sviluppo del metodo resta limi-
tata).
Teorema. Sia Φ(x, y(x);h; f) ∈ C0(D), D = S × [x0, x0 + β],
0 < h ≤ β, e lipschitziana in y. Allora un metodo one-step
esplicito e convergente se e solo se e consistente. Inoltre,
se il metodo e di ordine p, si ha
|ei| = |y(xi)− yi| ≤ k · hp
dove k e una costante indipendente da i e da h.
50
Propagazione degli errori di arrotondamento
Se indichiamo con ηi+1 l’errore di arrotondamento che si produce
nel calcolo di yi+1 ad ogni passo, si puo scrivere
yi+1 = yi + hΦ(xi, yi;h; f) + ηi+1
Se |ηi| ≤ η, ∀i, per l’errore globale di troncamento si ha la limitazione
|ei| ≤eL(xi−x0) − 1
L
(Cp hp +
η
h
)
In particolare, per il metodo di Eulero si ha
|ei| ≤eL(xi−x0) − 1
L
(M
2h+
η
h
)M = max
x∈[x0,x0+β]|y′′(x)|
51
Errore globale:
|ei| ≤eL(xi−x0) − 1
L
(M
2h+
η
h
)M = max
x∈[x0,x0+β]|y′′(x)|
hott
h
Errore globale
Errore di troncamento
Errore di arrotondamento
In corrispondenza del valore
ottimale hott =
√2η
M, dato
dall’intersezione della retta M2 h e
la retta ηh, l’errore di troncamento
e uguale a quello di arroton-
damento e la maggiorazione
dell’errore ha un minimo:
• per h < hott predomina l’errore
di arrotondamento;
• per h > hott predomina l’errore
di troncamento.
52
Esempio
Problema di Cauchy: f(x, y) = y(x)
{y′(x) = y(x) x ∈ [0,1]y(0) = 1
Soluzione esatta: y(x) = ex
Metodo di Eulero:
{yi+1 = yi + hf(xi, yi) = yi + hyi i = 1,2, . . . , ny0 = 1
1) x(1)i = ih1 i = 0,1, . . . , n1 h1 = 0.1 n1 = 10
2) x(2)i = ih2 i = 0,1, . . . , n2 h2 = 0.05 n2 = 20
53
0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 11
1.2
1.4
1.6
1.8
2
2.2
2.4
2.6
2.8
x
y y(x)
yi(2)
yi(1)
0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 10
0.02
0.04
0.06
0.08
0.1
0.12
0.14
x
e(1)i
=|y(xi)−y(1)
i|
e(2)i
=|y(xi)−y(2)
i|
Nota. M = maxx∈[0,1]
|y′′(x)| = maxx∈[0,1]
ex = e ' 2.7183
Se η = 0.5 · 10−14 ⇒ hott =
√2η
M' 6.06 · 10−8
54
Esercizio
(Primo esercizio della lezione)
{y′ = |y − x|, x ∈ [0,2]y(0) = 2
55
Soluzione
Ricordando che la soluzione esatta e y∗(x) = ex + x + 1, scegliendo il
passo h = 0.1 e risolvendo numericamente usando il Metodo di Eulero:
xi yi y∗i ei = |yi − y∗i |0.0 2.0000000000 2.00000000000000 0.000000000000000.1 2.2000000000 2.20517091807565 0.005170918075650.2 2.4100000000 2.42140275816017 0.011402758160170.3 2.6310000000 2.64985880757600 0.018858807576000.4 2.8641000000 2.89182469764127 0.027724697641270.5 3.1105100000 3.14872127070013 0.038211270700130.6 3.3715610000 3.42211880039051 0.050557800390510.7 3.6487171000 3.71375270747048 0.065035607470480.8 3.9435888100 4.02554092849247 0.081952118492470.9 4.2579476910 4.35960311115695 0.101655420156951.0 4.5937424601 4.71828182845905 0.12453936835904
56
Con h = 0.05
xi yi y∗i ei0.00 2.00000000000000 2.00000000000000 00.05 2.10000000000000 2.10127109637602 0.001271096376020.10 2.20250000000000 2.20517091807565 0.002670918075650.15 2.30762500000000 2.31183424272828 0.004209242728280.20 2.41550625000000 2.42140275816017 0.005896508160170.25 2.52628156250000 2.53402541668774 0.007743854187740.30 2.64009564062500 2.64985880757600 0.009763166951000.35 2.75710042265625 2.76906754859326 0.011967125937010.40 2.87745544378906 2.89182469764127 0.014369253852210.45 3.00132821597852 3.01831218549017 0.016983969511650.50 3.12889462677744 3.14872127070013 0.019826643922690.55 3.26033935811631 3.28325301786740 0.022913659751080.60 3.39585632602213 3.42211880039051 0.026262474368380.65 3.53564914232324 3.56554082901390 0.029891686690660.70 3.67993159943940 3.71375270747048 0.033821108031080.75 3.82892817941137 3.86700001661267 0.038071837201310.80 3.98287458838194 4.02554092849247 0.042666340110530.85 4.14201831780103 4.18964685192599 0.047628534124960.90 4.30661923369108 4.35960311115695 0.052983877465870.95 4.47695019537564 4.53570965931585 0.058759463940211.00 4.65329770514442 4.71828182845905 0.06498412331463
57
Riducendo il passo di discretizzazione, l’errore si e ridotto di circa la
meta .
Vogliamo ora stimare con quale valore del passo h l’errore di stima della
soluzione del PC nel punto x = 1 e sicuramente minore di 0.5 · 10−3.
Supponendo di trascurare gli errore di arrotondamento, per la stima
dell’errore globale di troncamento nel punto xi si puo usare il seguente
risultato
|ei| <M
2heL(xi−x0) − 1
L,
con M = maxx∈[0,1] |y′′(x)|, L la costante di Lipschitz della funzione
f(x, y).
58
Nel nostro caso, xi = 1, x0 = 0, L = 1 e y′′(x) = ex da cui M = e.
Ne risulta che
ei <e
2he− 1
1< 0.5 · 10−3
m
h <2 · 0.5 · 10−3
e(e− 1)=
10−3
e(e− 1)≈ 2.141 · 10−4 = 0.000214.
Infatti, scegliendo h = 0.0002, l’errore su y(1) e 0.2718 · 10−3.
59
Verificare questo risultato usando la funzione eulero.m
f=@(x,y)[abs(y-x)];
y_true = @(x)[exp(x)+x+1];
h = 0.0002;
n_step = 1/0.0002
n_step =
5000
[T]=eulero(0,2,f,h,n_step);
y_true_val = y_true(T(1,:));
ERR = abs(y_true_val-T(2,:));
ERR(end)
ans =
2.717783571872801e-004
60
Utilizziamo ora il metodo di Heun.
Essendo un metodo del secondo ordine, l’errore di troncamento globale
e
|ei| < C2h2 e
L(xi−x0) − 1
L,
con C2 opportuna costante.
61
function [T]=heun(x0,y0,fun,h,n_step)
%
% [T]= heun(x0,y0,fun,h,n_step)
% Soluzione di unequazione differenziale
% del primo ordine con il metodo di Heun
%
% Dati di input:
% x0, y0: condizione iniziale
% n: numero dei passi da eseguire
% h: passo di discretizzazione
% fun: espressione di f(x,y)
%
% OUTPUT
% T = matrice di dimensione 2xn_step
% la prima riga contiene il vettore dei nodi,
% la seconda riga contiene il vettore delle approssimazioni
% della soluzione nei nodi
%
62
xi(1) = x0;
yi(1) = y0;
for i = 2:n_step+1
fi = fun(xi(i-1),yi(i-1));
yi(i) = yi(i-1) + h/2*(fi+ fun(xi(i-1)+h,yi(i-1)+h*fi));
xi(i) = x0 + (i-1)*h;
end
T = [xi;yi];
63
Usando lo stesso passo h = 0.0002 si ha
Dal Command Window
f=@(x,y)[abs(y-x)];
y_true = @(x)[exp(x)+x+1];
h = 0.0002;
n_step = 1/h;
[T]=heun(0,2,f,h,n_step);
y_true_val = y_true(T(1,:));
ERR = abs(y_true_val-T(2,:));
ERR(end)
ans =
1.811916483518417e-008
L’errore in corrispondenza del punto x = 1 si e ridotto di ben 4 ordini
di grandezza.
64
In questo caso, invece, con
>> h = 0.025
h =
0.02500000000000
>> n_step = 1/h
n_step =
40
>> [T]=heun(0,2,f,h,n_step);
>> N=g(T(1,:));
>> E = abs(N-T(2,:));
>> E(end)
ans =
2.778840880699462e-004
si ha un errore confrontabile con quello ottenuto usando il metodo diEulero.
65
xi yi yi* ei = |yi*-yi|0 2.00000000000000 2.00000000000000 0
0.02500000000000 2.05031250000000 2.05031512052443 0.000002620524430.05000000000000 2.10126572265625 2.10127109637602 0.000005373719770.07500000000000 2.15287588626099 2.15288415088463 0.000008264623650.10000000000000 2.20515961963197 2.20517091807565 0.000011298443680.12500000000000 2.25813397250390 2.25814845306683 0.000014480562920.15000000000000 2.31181642618291 2.31183424272828 0.000017816545380.17500000000000 2.36622490447066 2.36624621661236 0.000021312141700.20000000000000 2.42137778486507 2.42140275816017 0.000024973295090.22500000000000 2.47729391004447 2.47732271619186 0.000028806147390.25000000000000 2.53399259964247 2.53402541668774 0.000032817045270.27500000000000 2.59149366232092 2.59153067486762 0.000037012546700.30000000000000 2.64981740814842 2.64985880757600 0.000041399427580.32500000000000 2.70898466129218 2.70903064598075 0.000045984688570.35000000000000 2.76901677303114 2.76906754859326 0.000050775562120.37500000000000 2.82993563509849 2.82999141461820 0.000055779519710.40000000000000 2.89176369336192 2.89182469764127 0.000061004279350.42500000000000 2.95452396185014 2.95459041966338 0.000066457813240.45000000000000 3.01824003713447 3.01831218549017 0.000072148355700.47500000000000 3.08293611307444 3.08301419748578 0.000078084411340.50000000000000 3.14863699593664 3.14872127070013 0.000084274763490.52500000000000 3.21536811989628 3.21545884837909 0.000090728482810.55000000000000 3.28315556293116 3.28325301786740 0.000097454936240.57500000000000 3.35202606311785 3.35213052691404 0.000104463796190.60000000000000 3.42200703534052 3.42211880039051 0.000111765049990.62500000000000 3.49312658842258 3.49324595743222 0.00011936900964
66
0.65000000000000 3.56541354269203 3.56554082901390 0.000127286321870.67500000000000 3.63889744799142 3.63903297596985 0.000135527978430.70000000000000 3.71360860214370 3.71375270747048 0.000144105326780.72500000000000 3.78957806988546 3.78973109996649 0.000153030081020.75000000000000 3.86683770227944 3.86700001661267 0.000162314333240.77500000000000 3.94542015661839 3.94559212718344 0.000171970565060.80000000000000 4.02535891683279 4.02554092849247 0.000182011659680.82500000000000 4.10668831441512 4.10688076532930 0.000192450914180.85000000000000 4.18944354987375 4.18964685192599 0.000203302052240.87500000000000 4.27366071472993 4.27387529396710 0.000214579237170.90000000000000 4.35937681407153 4.35960311115695 0.000226297085420.92500000000000 4.44662978967772 4.44686826035815 0.000238470680430.95000000000000 4.53545854372894 4.53570965931585 0.000251115586910.97500000000000 4.62590296311707 4.62616721098261 0.000264247865531.00000000000000 4.71800394437098 4.71828182845905 0.00027788408807
67
k1 = f(x0, y0) = |y0 − x0| = |2− 0| = 2
k2 = f(x0+h, y0+h·k1) = |y0+h·k1−x0−h| = |2+0.025·2−0.025| = 2.025
y1 = y0 +h
2(k1 + k2) = 2 + 0.025/2 · (2 + 2.025) = 2.0503125
68
x1 = x0 + h = 0.025
k1 = f(x1, y1) = |2.0503125− 0.025| = 2.0253125
k2 = f(x1 + h, y1 + h · k1) = |y1 + h · k1− x1 − h| =
= |2.0503125 + 0.025 · 2.0253125− 0.025− 0.025| = 2.0509453125
y2 = y1 +h
2(k1 + k2) =
= 2.0503125+0.025/2·(2.0253125+2.0509453125) = 2.10126572265625
69
Aumentando l’ordine del metodo numerico, per esempio scegliendo un
Runge-Kutta del 4o ordine
function [T]=RK4(x0,y0,fun,h,n_step)
%
% [T]= RK4(x0,y0,fun,h,n_step)
% Soluzione di unequazione differenziale
% del primo ordine con il metodo di Runge-Kutta di ordine 4
%
% INPUT:
% x0, y0: condizione iniziale
% n: numero dei passi da eseguire
% h: passo di discretizzazione
% fun: espressione di f(x,y)
% OUTPUT
% T = matrice di dimensione 2xn_step
% la prima riga contiene il vettore dei nodi,
% la seconda riga contiene il vettore delle approssimazioni
% della soluzione nei nodi
70
xi(1) = x0;
yi(1) = y0;
for i = 2:n_step+1
k1 = fun(xi(i-1),yi(i-1));
k2 = fun(xi(i-1)+h/2,yi(i-1)+h/2*k1);
k3 = fun(xi(i-1)+h/2,yi(i-1)+h/2*k2);
k4 = fun(xi(i-1)+h,yi(i-1)+h*k3);
yi(i) = yi(i-1) + h/6*(k1+2*k2+2*k3+k4);
xi(i) = x0 + (i-1)*h;
end
T = [xi;yi];
71
e usando h = 0.025, si ha
>> [T]=RK4(0,2,f,h,n_step);
>> N=g(T(1,:));
>> E = abs(N-T(2,:));
>> E(end)
ans =
8.666188655581664e-009
L’errore piu piccolo di quasi 5 ordini di grandezza.
72
Scegliendo h = 0.5 (si eseguono solo due passi del metodo!!!)
>> h = 0.5
h =
0.50000000000000
>> n_step = 1/h
n_step =
2
>> [T]=RK4(0,2,f,h,n_step);
>> N=g(T(1,:));
>> E = abs(N-T(2,:));
>> E(end)
ans =
9.356370527955349e-004
73
xi yi yi* ei = |yi*-yi|
0 2.00000000000000 2.00000000000000 00.50000000000000 3.14843750000000 3.14872127070013 0.000283770700131.00000000000000 4.71734619140625 4.71828182845905 0.00093563705280
k1 = f(x0, y0) = 2
k2 = f(x0 + h/2, y0 + h/2 · k1) = |2 + 0.5/2 · 2− 0.5/2| = 2.25
k3 = f(x0 + h/2, y0 + h/2 · k2) = |2 + 0.5/2 · 2.25− 0.5/2| = 2.3125
k4 = f(x0 + h, y0 + h · k3) = |2 + 0.5 · 2.3125− 0.5| = 2.65625
y1 = y0 + h/6 · (k1 + 2 · k2 + 2 · k3 + k4) =
= 2 + 0.5/6 · (2 + 2 · 2.25 + 2 · 2.3125 + 2.65625) = 3.1484375
74
x1 = x0 + h = 0.5
k1 = f(x1, y1) = |3.1484375− 0.5| = 2.6484375
k2 = f(x1+h/2, y1+h/2·k1) = |3.1484375+0.5/2·2.6484375−0.5−0.5/2| = 3.060546875
k3 = f(x1+h/2, y1+h/2·k2) = |3.1484375+0.5/2·3.060546875−0.75| = 3.16357421875
k4 = f(x1 +h, y1 +h ·k3) = |3.1484375 + 0.5 ·3.16357421875−1| = 3.730224609375
y2 = y1 + h/6 ∗ (k1 + 2 · k2 + 2 · k3 + k4) =
3.1484375+0.5/6·(2.6484375+2·3.060546875+2·3.16357421875+3.730224609375) =
= 4.71734619140625
x2 = x0 + 2 · h = 1
75
Esercizio
Si integri la seguente equazione
y′ = 3y − 4e−x, y(0) = 1
nell’intervallo [0,10], usando il metodo di Runge-Kutta del quarto
ordine con passo di discretizzazione h = 0.1.
Si confronti il risultato con la soluzione analitica y = e−x.
76
Usando la funzione RK4.m si ha
>> F =@(x,y)[3*y - 4*exp(-x)];
>> [T] = RK4(0,1,f,0.1,100);
>> T(:,1:20:end)’
0.0000e+000 1.0000e+000
2.0000e+000 1.3250e-001
4.0000e+000 -1.1237e+000
6.0000e+000 -4.6056e+002
8.0000e+000 -1.8575e+005
1.0000e+001 -7.4912e+007
Il comportamento della soluzione ottenuta e ”strano”. Infatti, guardando
alla soluzione analitica, ci si aspetterebbe di avere valori monotoni de-
crescenti delle yi, che, invece, prima decrescono e poi crescono in valore
assoluto (addirittura diventano negative!).
77
La soluzione generale dell’equazione differenziale e y(x) = Ce3x + e−x(si puo verificare per sostituzione).
La condizione iniziale y(0) = 1 implica C = 0 e quindi la soluzione e
proprio y = e−x.
La divergenza della soluzione prodotta e dovuta alla propagazione degli
errori di arrotondamento. Infatti, supponendo che la condizione iniziale
sia affetta da un piccolo errore ε, cioe y(0) = 1+ε, la soluzione analitica
dell’equazione differenziale diventa
y(x) = εe3x + e−x.
Quando x cresce, il termine dipendente da ε diventa dominante cau-
sando l’instabilita del metodo (sensibilita della soluzione alla condizione
iniziale).
78
Metodi impliciti
Un altro modo di derivare metodi numerici per la soluzione di equazionidifferenziali ordinarie e :
∫ xi+h
xiy′(x)dx
︸ ︷︷ ︸y(xi+h)−y(xi)
=∫ xi+h
xif(x, y(x))dx
e valutare l’integrale al secondo membro con una formula di quadraturausando un numero di nodi N con passo inferiore ad h,
⇒ y(xi + h)− y(xi) =N∑
r=1
crf(xi+h−lr, y(xi+h−lr)) +R(xi, y(xi);h; f)
⇒ yi+h e una combinazione lineare dei valori di f nei nodi e di yi
Esempio: usando la formula del trapezio
yi+1 = yi +h
2
(f(xi, yi) + f(xi+1, yi+1)
)Metodo di Crank-Nicolson
(metodo implicito del secondo ordine)
79
Metodi predictor-corrector
La soluzione viene approssimata con un metodo iterativo (corrector)in cui l’ approssimazione iniziale e ottenuta con un metodo esplicito(predictor). Per esempio,
Predictor : yi+1 = yi + hf(ti, yi) i ≥ 0
Corrector : y(r)i+1 = yi + h
2
(f(ti+1, y
(r−1)i+1 ) + f(ti, yi)
)1 ≤ r ≤ N
yi+1 = y(N)i+1
in cui e stato usato il metodo di Eulero esplicito come predittore e ilmetodo di Crank-Nicolson come correttore.
Oss:
1. La correzione si puo iterare N volte, con N scelto a priori oppure determinatousando un opportuno criterio di arresto.
2. Il metodo converge se h < 2L
, con L costante di Lipschitz di f .
3. Se L e grande (problemi stiff) la limitazione su h puo essere troppo restrittiva
4. In generale, l’ordine p di un metodo predictor-corrector e p = min{pP + 1, pC},con pP e pC gli ordini rispettivamente del predittore e del correttore.
80
Sistemi di equazioni differenziali
Sistema di equazioni differenziali del primo ordine
y′1(x) = f1(x, y1(x), . . . , yn(x))y′2(x) = f2(x, y1(x), . . . , yn(x)).................................................y′n(x) = fn(x, y1(x), . . . , yn(x))
y1(x0) = y10y2(x0) = y20.....................yn(x0) = yn0
Condizioni iniziali
In forma vettoriale:
{Y ′(x) = F (x, Y (x))Y (x0) = Y0
doveY (x) = [y1(x), y2(x), . . . , yn(x)]T
Y0 = [y1(x0), y2(x0), . . . , yn(x0)]T
F (x, Y (x)) = [f1(x, Y (x)), . . . , fn(x, Y (x))]T
Caso particolare: n = 2{y′(x) = f(x, y(x), z(x))z′(x) = g(x, y(x), z(x))
{y(x0) = y0z(x0) = z0
Condizioni iniziali
Nota. Per i sistemi di equazioni differenziali del primo ordine valgono teoremidi esistenza e unicita della soluzione analoghi a quelli validi nel caso diun’unica equazione differenziale.
81
Soluzione numerica disistemi di equazioni differenziali
{Y ′(x) = F (x, Y (x))Y (x0) = Y0
doveY (x) = [y1(x), y2(x), . . . , yn(x)]T
Y0 = [y1(x0), y2(x0), . . . , yn(x0)]T
F (x, Y (x)) = [f1(x, Y (x)), . . . , fn(x, Y (x))]T
Metodo di Eulero:
Yi+1 = Yi + hF (xi, Yi)(i = 0,1, . . . , n)
Y0 = Y (x0)dove
Y0(X0) = [y1(x0), y2(x0), . . . , yn(x0)]T
Yi = [y1i, y2i, . . . , yni]T
F (xi, Yi) = [f1(xi, Yi), . . . , fn(xi, Yi)]T
Caso particolare: n = 2
y′(x) = f(x, y(x), z(x))z′(x) = g(x, y(x), z(x))y(x0) = y0z(x0) = z0
→
yi+1 = yi + hf(xi, yi, zi) i = 0,1, . . . , nzi+1 = zi + hg(xi, yi, zi)y0 = y(x0)z0 = z(x0)
82
Metodo di Heun:Yi+1 = Yi +
h
2[F (xi, Yi) + F (xi + h, Yi + hF (xi, Yi))] i = 0,1, . . . , n
Y0 = Y (x0)
doveY0 = [y1(x0), y2(x0), . . . , yn(x0)]T
Yi = [y1i, y2i, . . . , yni]T
F (xi, Yi) = [f1(xi, Yi), . . . , fn(xi, Yi)]T
Caso particolare: n = 2
y′(x) = f(x, y(x), z(x))z′(x) = g(x, y(x), z(x))y(x0) = y0z(x0) = z0
yi+1 = yi +h
2[f(xi, yi, zi) + f(xi + h, yi + hf(xi, yi, zi)), zi + hg(xi, yi, zi))]
zi+1 = zi +h
2[g(xi, yi, zi) + g(xi + h, yi + hf(xi, yi, zi)), zi + hg(xi, yi, zi))]
i = 0,1, . . . , ny0 = y(x0)z0 = z(x0)
83
Metodo di Runge-Kutta classico (IV ordine):
Yi+1 = Yi +h
6[K1 + 2K2 + 2K3 +K4]
i = 0,1, . . . , nY0 = Y (x0)
dove
K1 = F (xi, Yi)
K2 = F
(xi +
h
2, Yi +
h
2K1
)
K3 = F
(xi +
h
2, Yi +
h
2K2
)
K4 = F (xi + h, Yi + hK3)
Caso particolare: n = 2
y′(x) = f(x, y(x), z(x))z′(x) = g(x, y(x), z(x))y(x0) = y0
z(x0) = z0
yi+1 = yi +h
6[k1 + 2k2 + 2k3 + k4] zi+1 = zi +
h
6[t1 + 2t2 + 2t3 + t4]
(i = 0,1, . . . , n)y0 = y(x0) z0 = z(x0)
k1 = f (xi, yi, zi) t1 = g (xi, yi, zi)
k2 = f
(xi +
h
2, yi +
h
2k1, zi +
h
2t1
)t2 = g
(xi +
h
2, yi +
h
2k1, zi +
h
2t1
)
k3 = f
(xi +
h
2, yi +
h
2k2, zi +
h
2t2
)t3 = g
(xi +
h
2, yi +
h
2k2, zi +
h
2t2
)
k4 = f (xi + h, yi + hk3, zi + ht3) t4 = g (xi + h, yi + hk3, zi + ht3)
84
Esercizio
Si consideri la seguente equazione differenziale in z(t)
z′′+ 2(z′)2 + 4z = t
con condizioni iniziali
z(0) = 0, z′(0) = 1
e lo si risolva numericamente usando due passi del metodo di Eulero
con passo di discretizzazione h = 1/2.
85
Definendo il vettore
y =
(y1y2
)=
(zz′)
l’equazione precedente si riduce ad un sistema di due equazioni dif-
ferenziali ordinarie del primo ordine
{y′1 = y2y′2 = t− 2y2
2 − 4y1
con condizioni iniziali
{y1(0) = 0y2(0) = 1
86
e che in forma vettoriale si scrive
dy
dt= f(y, t) y(0) =
(01
)
con
y =
(y1y2
)=
(zz′), e f =
(y2
t− 2y22 − 4y1
).
Usando il metodo di Eulero calcoliamo
yj+1 = yj + hf(tj,yj)
87
cioe
t0 = 0 y1 = y0+hf(t0,y0) =
(01
)+0.5·
(1
0− 2 · 12 − 4 · 0
)=
(1/2
0
)
t1 = 0.5 y2 = y1+hf(t1,y1) =
(1/2
0
)+0.5·
(0−3/2
)=
(1/2−3/4
)
88
function [T] = eulero_sistemi(x0,y0,z0,h,n_step,ffun,gfun)
% Soluzione di un sistema di due equazioni differenziali
% del primo ordine con il metodo di Eulero
%
% INPUT
% x0, y0,z0 = consizioni iniziali
% h = passo di discretizzazione
% n = numero passi
% ffun = input(funzione f(x,y,z): );
% gfun = input(funzione g(x,y,z): );
%
% OUTPUT
% T = matrice contenente [xi;yi;zi]
89
xi(1) = x0; yi(1) = y0; zi(1) = z0;
for i = 2:n_step+1
x = xi(i-1);
y = yi(i-1);
z = zi(i-1);
k1 = ffun(x,y,z);
t1 = gfun(x,y,z);
xi(i) = x + h;
yi(i) = y + h*k1;
zi(i) = z + h*t1;
end
figure(1)
plot(xi,yi,’r’,xi,zi,’b’)
T = [xi;yi;zi];
90
Dal Command window
>> ffun = @(t,y,z)[z];
>> gfun = @(t,y,z)[t-2*z^2-4*y];
>> [T] = eulero_sistemi(0,0,1,.5,2,ffun,gfun);
>> disp([’ ti ’,’ y1(i) ’, ’ y2(i) ’ ]), disp(T’)
ti y1(i) y2(i)
0 0 1.0000
0.5000 0.5000 0
1.0000 0.5000 -0.7500
Usando il metodo di RK del quarto ordine con lo stesso passo di di-scretizzazione si ha
>> [T] = RK_sistemi(0,0,1,.5,2,ffun,gfun);
>> disp([’ ti ’,’ y1(i) ’, ’ y2(i) ’ ]), disp(T’)
ti y1(i) y2(i)
0 0 1.0000
0.5000 0.3148 0.3259
1.0000 0.3761 -0.0401
91
Usando il metodo di Eulero con passo di discretizzazione h = 0.05:
>> [T] = eulero_sistemi(0,0,1,.05,20,ffun,gfun);>> disp([’ ti ’,’ y1(i) ’, ’ y2(i) ’ ]), disp(T’)
ti y1(i) y2(i)0 0 1.0000
0.0500 0.0500 0.90000.1000 0.0950 0.81150.1500 0.1356 0.73160.2000 0.1722 0.65850.2500 0.2051 0.59070.3000 0.2346 0.52730.3500 0.2610 0.46760.4000 0.2844 0.41100.4500 0.3049 0.35720.5000 0.3228 0.30600.5500 0.3381 0.25710.6000 0.3509 0.21040.6500 0.3614 0.16570.7000 0.3697 0.12320.7500 0.3759 0.08270.8000 0.3800 0.04440.8500 0.3823 0.00820.9000 0.3827 -0.02580.9500 0.3814 -0.05741.0000 0.3785 -0.0865
92
Mentre il metodo di RK con passo h = 0.05 restituisce
>> [T] = RK_sistemi(0,0,1,.05,20,ffun,gfun);>> disp([’ ti ’,’ y1(i) ’, ’ y2(i) ’ ]), disp(T’)
ti y1(i) y2(i)0 0 1.0000
0.0500 0.0476 0.90570.1000 0.0907 0.82110.1500 0.1298 0.74400.2000 0.1652 0.67300.2500 0.1972 0.60700.3000 0.2260 0.54520.3500 0.2518 0.48700.4000 0.2748 0.43190.4500 0.2950 0.37970.5000 0.3128 0.33000.5500 0.3281 0.28270.6000 0.3411 0.23780.6500 0.3519 0.19500.7000 0.3606 0.15440.7500 0.3674 0.11600.8000 0.3723 0.07980.8500 0.3754 0.04580.9000 0.3769 0.01420.9500 0.3768 -0.01491.0000 0.3754 -0.0415
93
Esercizio
Consideriamo l’oscillatore armonico descritto dall’equazione differen-
ziale del secondo ordine (n = 2):
u′′(x) + 2αu′(x) + β2u(x) = 0
↓u′′(x) = −2αu′(x)− β2u(x) := g(x, u(x), u′(x))
u′′(x) = g(x, u(x)︸ ︷︷ ︸y(x)
, u′(x)︸ ︷︷ ︸z(x)
)
u(x0) = u0u′(x0) = u1
→
y(x) = u(x)
y′(x) = z(x)z′(x) = g(x, y(x), z(x))y(x0) = y0z(x0) = y1
Per risolvere il sistema differenziale ottenuto, si possono applicare i
metodi gia visti per la soluzione dei sistemi (ad esempio, il metodo di
Runge-Kutta classico).
94
Script MATLAB% Soluzione di un sistema di due equazioni differenziali% del primo ordine con il metodo di Runge-Kutta classico%% Input%clearx0 = input(’x0: ’);y0 = input(’y0: ’);z0 = input(’z0: ’);h = input(’passo: ’);n = input(’numero passi: ’);ffun = input(’funzione f(x,y,z): ’);f = inline(ffun,’x’,’y’,’z’);gfun = input(’funzione g(x,y,z): ’);g = inline(gfun,’x’,’y’,’z’);%% Algoritmo%xi(1) = x0; yi(1) = y0; zi(1) = z0;for i = 2:n
x = xi(i-1);y = yi(i-1);z = zi(i-1);k1 = f(x,y,z);t1 = g(x,y,z);
95
k2 = f(x+0.5*h,y+0.5*h*k1,z+0.5*h*t1);t2 = g(x+0.5*h,y+0.5*h*k1,z+0.5*h*t1);k3 = f(x+0.5*h,y+0.5*h*k2,z+0.5*h*t2);t3 = g(x+0.5*h,y+0.5*h*k2,z+0.5*h*t2);k4 = f(x+h,y+h*k3,z+h*t3);t4 = g(x+h,y+h*k3,z+h*t3);xi(i) = x + h;yi(i) = y + h*(k1+2*k2+2*k3+k4)/6;zi(i) = z + h*(t1+2*t2+2*t3+t4)/6;
endfigure(1)plot(xi,yi,’r’,xi,zi,’b’)
0 5 10 15 20 25 30−1
−0.8
−0.6
−0.4
−0.2
0
0.2
0.4
0.6
0.8
1
x
y
Oscillazioni smorzate: α=0.1, β=1
y(x)=u(x)
z(x)=u’(x)
0 1 2 3 4 5 6 7 8 9 10−0.4
−0.2
0
0.2
0.4
0.6
0.8
1
x
y
Smorzamento supercritico: α=1, β=1
y(x)=u(x)
z(x)=u’(x)
Dati del problema Dati del problema
x0: 0 x0: 0y0: 1 y0: 1z0: 0 z0: 0passo: 0.2 passo: 0.2numero passi: 150 numero passi: 150funzione f(x,y,z): ’z’ funzione f(x,y,z): ’z’funzione g(x,y,z): ’-2*0.1*z-y’ funzione g(x,y,z): ’-2*z-y’
96
Esercizio
1. Utilizzare il programma per approssimare la soluzione del problema
del pendolo non lineare. Confrontare graficamente la soluzione
approssimata ottenuta con quella del problema lineare.
2. Risolvere il sistema di equazioni differenziali del modello preda-
predatore, per i seguenti valori dei parametri k1 = k2 = 1, µ1 =
300, µ2 = 200, y10 = 400, y20 = 100, usando il metodo di Eulero
con passi 0.1, 0.01 e 0.001, il metodi di Runge Kutta del quarto
ordine con passi 0.1 e 0.01 e il metodo di Crank-Nicolson con passi
0.01 e 0.001. Confrontare i risultati
97
Riferimenti bibliografici
L. Gori, Calcolo Numerico: Cap. 9 §§ 9.1-9.6, 9.14 (esclusi metodi di Adams-Bashfort, Milne e Numerov)
L. Gori, M.L. Lo Cascio, F. Pitolli, Esercizi di Calcolo Numerico: 6.1, 6.2, 6.3, 6.4,
6.5, 6.7, 7.76, 7.80, 7.85
98
Esercizi
99
Esercizio 1
Dato il problema di Cauchy{
y′ = 1 + y2, x ∈ [0,1]y(0) = 0
dopo aver verificato che y = tan(x) e la soluzione esatta,
1. risolvere numericamente usando il metodo di Eulero esplicito con
passo h = 12 e produrre una stima per y(1)
2. scegliendo un passo h = 12p e trascurando l’influenza degli errori di
arrotondamento, stabilire per quale valore di p l’errore globale di
troncamento per y(1) diventa minore di 0.5 · 10−4.
100
Soluzione
Si osserva che se y∗ = tan(x), allora (y∗)′ = 1cos(x)2.
Usando la relazione fondamentale della trigonometria si ha
(y∗)′ = cos2(x) + sin2(x)
cos2(x)= 1 + tan2(x) = 1 + (y∗)2.
Quindi, y∗ e soluzione del problema di Cauchy.
La soluzione e unica in [0,1] × R: Infatti, poiche risulta |fy| = |2y| =
|2 tanx| ≤ 2 tan(1), la funzione f(x, y) = 1+y2 e Lipschitziana rispetto
a y, uniformememnte rispetto a x nel dominio considerato.
101
Il metodo di Eulero produce approssimazioni della soluzione nei punti
xi = x0 + ih.
Se h = 1/2 e x0 = 0, la stima per y(1) si ottiene dopo due iterazioni:
x1 = 1/2 y(x1) = y(0) + h(1 + y2(0)) = 1/2
x2 = 1 y(x2) = y(x1)+1
2(1+y2(x1)) =
1
2+
1
2
(1 +
1
4
)=
9
8= 1.1250
Poiche y∗(1) = 1.5574, l’errore relativo e |y∗(1)−y(1)|y∗(1) = 1.5574−1.1250
1.5574 =
0.2776.
102
Una stima per l’errore globale di troncamento ei per il metodo di Euleroe data
|ei| <M
2heL(xi−x0) − 1
L,
con M = maxx∈[0,1] |y′′(x)|, L la costante di Lipschitz della funzionef(x, y) = 1 + y2 rispetto a y, xi il punto in cui valutare l’errore e h ilpasso adottato.
In questo caso h = 12p, xi − x0 = 1− 0 = 1,
M = maxx∈[0,1]
|2yy′| = maxx∈[0,1]
|2y(x)(1 + y(x)2)| =
= |2 tan(1)(1 + tan2(1))| = 10.6699
(e stata usata la monotonia di tan(x) nell’intervallo [0,1])
mentre come costante di Lipschitz si assume
L = max(x,y)
|fy(x, y)| = max(x,y)
|2y| = |2y(1)| = 3.1148.
103
Quindi,
p > log2
(10.6699
2
e3.1148 − 1
3.1148
1
0.5 · 10−4
)= 19.49,
ovvero p = 20, e quindi h = 1220 = 9.536743164062500 · 10−7.
Se non si trascurano gli errori di arrontondamento, l’ errore di arroton-damento si maggiora nel modo seguente
|ei| <(M
2h+
η
h
)eL(xi−x0) − 1
L,
dove η e tale che |ηi| ≤ |η|, ∀ i.
Allora e possibile definire il passo ottimo hott che rende minimo l’erroretotale (troncamento + propagazione), cioe
hott =
√2η
M
Supponendo che η = 0.5 · 10−15, si ha
hott =
√10−15
10.6699= 9.680991201464672 · 10−9
104
Esercizio 2
Si vuole risolvere il problema di Cauchy{y′(x) = y + cos(x) , x ∈ (−2,0) ,y(−2) = 0 ,
con il metodo di Eulero.
2.1) Sapendo che y(x) ∈ C(3)([−2,0]) e che |y(p)(x)| ≤ 4, p = 0,1,2,3,
per x ∈ [−2,0], determinare quale passo di integrazione garantisce
che l’errore globale di troncamento in x = 0 sia inferiore a 0.5·10−5.
2.2) Sapendo che ad ogni passo si ha un’errore di arrotondamento sulla
5 cifra decimale, individuare il passo ottimo.
105
Soluzione
2.1) Bisogna individuare un passo h tale che
|ei| ≤h
2Me(0−x0)L − 1
L≤ 0.5 · 10−5 ,
dove M ≥ |y′′(x)| (e quindi M = 4) e L e la costante di Lipschitzdi f rispetto a y. Poiche f(x, y) = y + cosx e
∣∣∣∣∣∂f(x, y)
∂y
∣∣∣∣∣ = |1| ⇒ L = 1,
si ha
h
24e2 − 1
1≤ 0.5 · 10−5 ,
da cui si ottiene
h ≤ 0.5 · 10−5
2(e2 − 1)≈ 3.9129 · 10−7.
Scegliendo h = 1 · 10−7, valore di poco sopra alla maggiorazioneottenuta, sono necessari 2 · 107 passi del metodo di Eulero perarrivare all’estremo destro dell’intervallo di integrazione.
106
2.2) Il passo ottimo e dato da
hott =
√2η
M=
√2 · 0.5 · 10−5
4≈ 0.0016 .
Il passo individuato al punto precedente, essendo molto piu piccolo
di hott, puo produrre un errore di propagazione elevato.
ESERCIZIO 3
Si consideri il problema di Cauchy{y′ = y − 2sin(x), 0 ≤ x ≤ 5y(0) = 1
.
3.1) Verificare che la soluzione esiste ed e unica in D = [0,5] ×R e che la funzioney(x) = cos(x) + sin(x) e la soluzione esatta.
3.2) Supponendo di voler risolvere numericamente il problema di Cauchy con il metododi Eulero,
– posto il passo di discretizzazione h = 52m e trascurando gli errori di arroton-
damento, stabilire per quale valore di m ∈ N l’errore globale di troncamentoper y(5) diventa minore di 0.5 · 10−3;
– supponendo di introdurre ad ogni passo un errore di arrontondamento infe-riore a 0.5 · 10−3, calcolare il passo di discretizzazione ottimo e confrontarlocon quello ottenuto al passo precedente.
107
ESERCIZIO 4
Si consideri il seguente problema di Cauchy (PC):{y′(t) = y(t)− t2 + 1 t ∈ [0,2]y(0) = 1/2
4.1 Scelto un passo h = 13m, con m ∈ N e sapendo che max0≤t≤2 |y′′(t)| ≤ 1.7, stabilire
per quali valori di m l’errore globale massimo della soluzione numerica del (PC)calcolata con il metodo di Eulero, trascurando gli errori di arrotondamento,risulta non superiore a 0.5 · 10−3.
4.2 Confrontare il passo h piu grande calcolato al punto precedente con il valore delpasso ottimo hott calcolato quando si considera un errore massimo di arrotonda-mento pari a 0.5 · 10−9.
108
ESERCIZIO 5
5.1 Illustrare il problema della soluzione numerica di una equazione dif-
ferenziale con particolare riferimento al metodo di Eulero esplicito.
5.2 Stabilire se il metodo di Eulero e adatto ad approssimare la soluzione
z(y) =(
23y)3
2, per y ≥ 0, del seguente problema di Cauchy
{z′ = z
13
z(0) = 0.
109
