Analisi Numerica

(A.A. 2014-2015)

Appunti delle lezioni:

Equazioni differenziali ordinarie

Docente Vittoria Bruni

Email: vittoria.bruni@sbai.uniroma1.it

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