1

LUISSCorso di

Metodi matematici per economia e finanza. Prof. F. Gozzi

Lezione del 16/11/2009Integrazione numerica di equazioni differenziali

- Prof.ssa G. Rotundo

2

Testi di riferimento

MATLAB – manuale di riferimento

I. Capuzzo Dolcetta, M. Falcone, L’analisi al calcolatore, Zanichelli, ISBN 88 – 08 – 03904 - 8

J. Stoer, Introduzione all’analisi numerica, Zanichelli ed., 1974

3

Il problema di Cauchy

inizialecondizione)(

evoluzione di legge))(,()('

00 Xxtx

ttxtftx T

Sotto opportune ipotesi è possibile dimostrare il teorema di esistenza ed unicità delle soluzioni.

N.B.: il teorema di esistenza ed unicità fornisce un risultato soltanto in merito ad esistenza ed unicità, ma non illustra alcun metodo per trovare la soluzione.

Domanda: come trovare le soluzioni?

Anche quando la soluzione esiste ed è unica può essere molto complesso, se non impossibile, determinare la sua espressione analitica. Diventa quindi estremamente importante conoscere alcuni metodi numerici per approssimare la soluzione.

4

Metodi numerici di integrazione

Basati sulla formula di Taylor

A. Metodo di Eulero 1. Analisi della qualità del metodo 2. Stima dell’errore globale di approssimazione a. Illustrazione tramite un esempio b. Stima (con dimostrazione) B. Alri metodi di ordine superiore: 1.Metodo di Eulero modificato 2.Metodo di Runge-Kutta

5

Il metodo di Eulero

Leonhard Euler (1707-1783)

6

Il metodo di Eulero

E’ il metodo più semplice.

Per illustrarlo si parte dalla definizione di derivata:

)(')()(

lim0

txh

txhtxh

Per h abbastanza piccolo

))(,()(')()(

txtftxh

txhtx

))(,()()( txtfhtxhtx Da cui

7

))(,()()( txtfhtxhtx

))(,()()( txtfhtxhtx

Questa relazione suggerisce di definire una soluzione

approssimata a partire dalla successione definita per ricorrenza da

Idea: fisso h e valuto l’espressione

sull’insieme di punti

hthtt

hthtt

htt

t

3

2

023

012

01

0

La quantità h viene chiamata lunghezza del passo oppure passo dello schema di approssimazione

8

In corrispondenza a tk si può definire una successione di xk :

1

2

1

0

kt

t

t

t

))(,()()( txtfhtxhtx

),()()(

),()()(

),()()(

)(

11

111122

000011

00

kkkkkk xtfhxhtxtxx

xtfhxhtxtxx

xtfhxhtxtxx

txx

9

Interpretazione del metodo di Eulero: il punto successivo è ottenuto dal precedente spostandosi con pendenza fissata per un intervallo di ampiezza h

10

Esempio

Considero la seguente funzione:

Si verifica che soddisfa

una equazione differenziale:

ttxtx

txtxt

ttxt

ttx

ttx

/)(2)('

)(2)('

2)('

2)('

)(

2

2

11

Esempio (sarà ripreso nella lezione in aula informatica)

Considero l’equazione differenziale insieme alle condizioni iniziali x(0.5)=0.25

0.5 1 1.5 2 2.50

1

2

3

4

5

6

Problema: la successione approssimante (spezzata in nero con i punti evidenziati) si allontana ben presto dalla soluzione .

E’ necessario dare una stima della bontà della approssimazione:

Ricordando che:

))(,()()( txtfhtxhtx x’(t)=f(t,x(t))

12

Sviluppo di Taylor arrestato al I ordine

Osservo che le due condizioni

sono descritte per t=0 dallo sviluppo di Taylor:

22

2

)(")0(')0()( Ch

sxhxhxhx

))(,()()( txtfhtxhtx

x’(t)=f(t,x(t))

),0(,2

)(")0(')0()( 2 hs

sxhxhxhx

2

)("sup

),0(

sxC

hs

13

Osservazione

Il metodo di Eulero viene detto metodo ad un passo perché il calcolo del punto successivo dipende unicamente dal punto precedente.

Un generico metodo ad un passo è definito da una successione per ricorrenza del tipo

),(

)0(

1

0

kkk xhhxx

xx

),(),( kk xhfxh Il metodo di Eulero è un caso particolare:

14

L’errore di approssimazione è piccolo se h è piccolo?

Analisi della qualità del metodo

Voglio determinare

1. se l’errore di approssimazione è piccolo se h è piccolo (quindi definisco una misura per l’errore e dimostro che il limite è 0 per h0)

2. l’errore globale di approssimazione (cosa succede se itero il procedimento N volte ed N∞)

15

1. Misura della qualità di un metodo ad un passo

Errore locale di discretizzazione in x(0) :

),(),(),( 000 hxhxhxe

0se)),0(,0(

0se,)0()(

),( 0

hxf

hh

xhxhx

16

Dimostro che 0)),0((lim0

hxeh

Svolgimento:

0lim

))(,()0()(lim

))(,()0()(lim

))(,()0()(

lim)),0((lim

2

00

0

00

h

Ch

h

txthfxhx

h

txthfxhx

txtfh

xhxhxe

hh

h

hh

17

Osservazioni

Poiché l’errore locale di discretizzazione è dello stesso ordine di h, per h 0, si dice che il metodo di Eulero è un metodo del primo ordine.

La dimostrazione della convergenza non garantisce che la soluzione approssimata sia buona, vediamo un ulteriore esempio.

18

2. Stima dell’errore globale di discretizzazione

a. Illustrazione del problema tramite un esempiob.stima

19

2.a Illustrazione del problema tramite un esempio

Considero la seguente funzione:

Si verifica che soddisfa

una problema di Cauchy: )()('

)('

)(

txtx

etx

etxt

t

1)0(

)()('

x

txtx

20

Sviluppo in serie di Taylor

1)0(

))(()()('

x

txftxtx

hxhxhxx 1)0(')0()(1

212 )1()1()1()1( hhhhhfhxx

Calcolo i punti successivi

Osservo che la variabile t non appare da sola, ma soltanto tramite x

Osservo che la funzione f(.) è la funzione identica (restituisce l’argomento, immutato).

322223 )1()1()1(1 hhhhhfhxx

21

Quindi la funzione è:tetx )(

khx k )1( e le iterate sono:

Voglio stimare l’errore (t=hk, così lo stimo nei punti della successione tk)

)1log()1()( hkkhkkhk eehexkhx

22

Promemoria: Teorema di Lagrange

))((')()( abcfafbf

Sia f continua in [a,b], derivabile in (a,b). Allora esiste un punto c in (a,b) tale che

23

Usando il teorema di Lagrange

0))1log((

))1log(()1log(

hhk

hkkheee hkkh

)),1(log( hh

k

kxkhx )(

kpk

x

xkhxp

)(

Con a=klog(1+h) e b=kh

Perché >0, quindi e >1

Si può dimostrare anche che:

24

2.b Errore globale di discretizzazione

Misura della qualità della approssimazione a tempi molto maggiori di quello in cui è data la condizione iniziale

kTt

xtxhE

)(sup)(],0[

25

Stima di E(h) per il metodo di Eulero

Ipotesi ulteriori: f ed f ‘ limitate in R ;)('sup;)(sup tfLtfMRtRt

Fissato T e scelto h=T/N, N intero positivo fissato, vogliamo dimostrare che

2/)1()( LTN eMhTxx

Osservazione 1: Questa stima dimostra che si può approssimare la soluzione esatta con la precisione desiderata su un intervallo [0,T] di ampiezza qualsiasi.

Osservazione 2: con le ipotesi fatte la soluzione esiste in [0,+∞)

Osservazione 3: la stima di questa disuguaglianza permette di dimostrare la convergenza a zero di |E(h)| per h0

Osservazione 4: se T è grande, h dovrà essere scelto molto piccolo se si vuole avere una buona approssimazione.

Tesi

26

• Svolto fino a qui il 16 nov 09

27

;)('sup;)(sup tfLtfMRtRt

2/)1()( LTN eMhTxx

dhkxhkfhkxkhxh

)))1((,)1(())1(()(0

Dimostrazione

Osservazione: se x è soluzione del PC-EDO, allora

stt

s

tsdsxsfsxdxfsxtx0

,,))(,()())(,()()(

)(khxxd kk Ponendo si ha dunque

dhkxhkfhkx

yhfydh

kkk

)))1((,)1(())1((

)(

0

11

Caso particolare:

Ometto la dipendenza da t

28

h

kk

h

kkk

dhkxfxfd

dhkxhkfhkx

yhfyd

0

11

0

11

)))1((()((

)))1((,)1(())1((

)(

Applicando il teorema di Lagrange e le ipotesi di limitatezza di f ed f ‘ si ha

h

kk

h

kk

h

kkk

dxhkxhkxhkxLd

dxhkxhkxhkxLd

dxhkxLdd

0

11

0

11

0

11

))1(())1(())1(((

))1(())1(())1(((

))1(((

|dk-1|Vd pagina successiva

29

Osservazione: per t=(k-1)h+ e s=(k-1)h l’equazione

stt

s

tsdsxsfsxdxfsxtx0

,,))(,()())(,()()(

diventa

Mdhkxfhkxhkx 0

)))1((())1(())1((

Pertanto

h

kkk dLMdLhdd0

11

kh

LMdLhd kk ,2

)1(2

1

Iterando questa relazione ‘all’indietro’ fino a k=0 e ricordando che d0=x(0)-y0=0 si trova:

Ipotesi di limitatezza

30

Iterando questa relazione ‘all’indietro’ fino a k=0 e

ricordando che d0=x(0)-y0=0 si trova:

2/]1)1[(

2/)1()1()1(1 122

k

kk

LhMh

LhLhLhLMhd

Scegliendo ora k=n e ricordando che h=T/N, osservando che (1+LT/N)N<eLT

2/)1()( LTN eMhTxx

31

Metodi di ordine superiore basati sullo sviluppo di Taylor

Vogliamo ora studiare dei metodi di approssimazione che, a parità di pazzo h, diano luogo a errori locali e globali di approssimazione più piccoli rispetto a quelli ottenuti per il metodo di Eulero.

In particolare, ci occuperemo di vari metodi ad un passo di ordine p, intendendo con questo che

ch

hte

hte

ph

h

),(lim

0),(lim

0

0 Ovvero che l’errore locale tende a zero con la stessa rapidità di hp per h0

In tal caso si dice che l’errore locale è un infinitesimo di ordine pe(t,h)=O(hp)

32

Utilizzando lo sviluppo di Taylor è possibile costruire facilmente dei metodi ad un passo di ordine comunque elevato, a patto che f sia derivabile un numero sufficiente di volte. Supponiamo dunque che f sia dotata di (p-1) derivate continue. Il suo sviluppo di Taylor in un intorno di t è

)(!/)(2/)('')(')()( )(2 hRptxhtxhthxtxhtx ppp

Da cui si ricava

hhRptxhthxtxh

txhtxp

pp /)(!/)(2/)('')(')()( )(1

33

Poiché x è soluzione del PC-EDO

))(()))(('())(())(('')('''

))(())((')('))((')(''

))(()('

22 txftxftxftxftx

txftxftxtxftx

txftx

6/)]())('()()(''[2/)()(')(),( 222 xfxfxfxfhxfxhfxfhx

E così via. Possiamo cioè esprimere tutte le derivate di y in funzione di f e delle sue derivate calcolate nel punto y(t).

Di conseguenza per avere un metodo di ordine 3, ad esempio, basterà definire

),(1 hxhxx kkk

E la successione come

34

Calcolo l’errore locale di discretizzazione in x(0) :

)(!3

)0("'

2

)0(")0(')0()( 3

32 hRx

hx

hxhxhx

),(),(),( 000 hxhxhxe

0se)),0(,0(

0se,)0()(

),( 0

hxf

hh

xhxhx

Ricordo anche l’espansione in serie di Taylor

Da cui

h

hRxh

xhxhxhx

)(

!3

)0("'

2

)0(")0(')0()( 332

Sostituendo l’espressione di (xo,h) :

),()(

!3

)0("'

2

)0(")0('),( 32 ht

h

hRxh

xhxhte

h

hRhxe

)(),( 3

0 )(),( 30 hOhxe

quindi

35

Osservazioni

1. Questi metodi possono essere usati per mdefinire metodi di approssimazione ad un passo di ordine comun que elevato.

2. I metodi ottenuti in questo modo richiedono la conoscenza di una formula esplicita di ciascuna delle derivate di f che compaiono in .

3. Siccome le derivate possono non essere facili, in alcuni casi si possono sostituire con i corrispondenti rapporti incrementali, ma questo appesantisce il calcolo.

36

Altri metodi di ordine superiore

Obbiettivo:

Ottenere metodi che abbiano ordine di convergenza elevato senza complicare troppo la trattazione.

Descriviamo la costruzione di alcuni metodi ad un passo di ordine 2 aventi queste caratteristiche.

37

Ipotesi: derivabile due volte

)()0,()0,(),( 2000 hOxhxhx h

)(!3

)0("'

2

)0(")0(')0()( 3

32 hRx

hx

hxhxhx

)()0,(')0,(

)(2

)()(')(),0(

200

2000

hOxhx

hOxfxf

hxfhxe

n

Per ottenere un metodo di ordine due basterà allora scegliere in modo tale che per ogni x si abbia

=0

38

Esempio

))(()(),( 000 xChfxBfxAfhx

0)0,(')0,(2

)()(')( 00

000 xhx

xfxfhxf n

)()(')()(2

)(')( 00

00 xfxhBCfxfBA

xfhxf

Con A, B, C parametri non negativi da determinare. La condizione

diventa

Che è verificata se A, B, C soddisfano il sistema di equazioni

2

11

BC

BA

39

Possibili scelte per A, B, C

Evidentemente le equazioni non determinano univocamente i parametri. A differenti soluzioni corrispondono differenti metodi ad un passo di ordine due.

Metodo di Heun: A=1/2, B=1/2, C=1

Metodo di Eulero modificato: A=0, B=1, C=1/2

40

Metodo di Runge-Kutta

Si ottiene definendo

6/)22(),0( 4321 kkkkhx

Con

)(

)2/(

)2/(

)(

304

203

102

01

hkxfk

hkxfk

hkxfk

xfk

Si può dimostrare che è un metodo ad un passo di ordine 4.

41

Osservazioni

1. Questi metodi hanno vantaggi rispetto a quelli basati sulla formula di Taylor. Per esempio si chiede solo che f abbia la derivata prima continua e sarà sufficiente calcolare f in un minor numero di punti.

2. Il risparmio del calcolo di f anche in un solo punto puo’ portare ad un grande risparmio nel tempo di calcolo complessivo in applicazioni che necessitano di un elevato numero di iterazioni

3. Sotto l’ipotesi che f sia dotata di derivate continue e limitate si puo’ dimostrare che E(h)=O(h2) per i metodi di Heun ed Eulero modificatoE(h)=O(h4) per il metodo di Runge-Kutta

42

Nota sui sistemi del primo ordine

I metodi visti finora si estendono ad un sistema di N equazioni del primo ordine.

NN

NN

xx

xx

xx

txftx

txftx

txftx

)0(

)0(

)0(

))(()('

))(()('

))(()('

22

11

22

11

Che in forma vettoriale diventa

0)0(

))(()('

xx

txftx

Lo spazio euclideo RN dove è assegnata la condizione iniziale e dove evolve la traiettoria x(t) prende il nome di spazio delle fasi