LUISS
Corso di
Metodi matematici per economia e finanza.
Prof. F. Gozzi
Lezione del 16/11/2009
Integrazione numerica di equazioni differenziali
- Prof.ssa G. Rotundo
1
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
2
Il problema di Cauchy
 x' (t )  f (t , x(t )) t  T

 x(t0 )  x0  X
legge di evoluzione
condizione iniziale
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 3
alcuni metodi numerici per approssimare la soluzione.
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
4
Il metodo di Eulero
Leonhard Euler (1707-1783)
5
Il metodo di Eulero
E’ il metodo più semplice.
Per illustrarlo si parte dalla definizione di derivata:
x(t  h)  x(t )
lim
 x' (t )
h 0
h
Per h abbastanza piccolo
x(t  h)  x(t )
 x' (t )  f (t , x(t ))
h
Da cui
x(t  h)  x(t )  h  f (t , x(t ))
6
x(t  h)  x(t )  h  f (t , x(t ))
Questa relazione suggerisce di definire una soluzione
approssimata a partire dalla successione definita per ricorrenza da
x(t  h)  x(t )  h  f (t , x(t ))
Idea: fisso h e valuto l’espressione
sull’insieme di punti
La quantità h viene chiamata
lunghezza del passo oppure passo
dello schema di approssimazione
t0
t1  t0  h
t 2  t1  h  t0  2h
t3  t 2  h  t0  3h

7
x(t  h)  x(t )  h  f (t , x(t ))
In corrispondenza a tk si può definire una successione di xk :
t0
x0  x(t0 )
t1 x1  x(t1 )  x(t0  h)  x0  h  f (t0 , x0 )
t 2 x2  x(t2 )  x(t1  h)  x1  h  f (t1 , x1 )


t k 1 xk 1  x(tk 1 )  x(tk  h)  xk  h  f (tk , xk )
8
Interpretazione del metodo di Eulero: il punto successivo è
ottenuto dal precedente spostandosi con pendenza fissata per un
intervallo di ampiezza h
9
Esempio
Considero la seguente
funzione:
Si verifica che soddisfa
una equazione
differenziale:
x(t )  t
x' (t )  2t
2
t  x' (t )  2t
t  x' (t )  2 x(t )
x' (t )  2 x(t ) / t
2
10
Esempio (sarà ripreso nella lezione
in aula informatica)
6
Considero l’equazione
differenziale insieme alle
condizioni iniziali
x(0.5)=0.25
Problema: la successione
approssimante (spezzata in nero con
i punti evidenziati) si allontana ben
presto dalla soluzione .
5
4
3
2
1
0
0.5
E’ necessario dare una stima della
bontà della approssimazione:
Ricordando che:
1
1.5
2
2.5
x(t  h)  x(t )  h  f (t , x(t ))
x’(t)=f(t,x(t))
11
Sviluppo di Taylor arrestato al I ordine
Osservo che le due condizioni
x(t  h)  x(t )  h  f (t , x(t ))
x’(t)=f(t,x(t))
sono descritte per t=0 dallo sviluppo di Taylor:
x" ( s )
x(h)  x(0)  h  x' (0)  h 
, s  (0, h)
2
2
x" ( s)
2
x(h)  x(0)  h  x' (0)   h 
 Ch
2
2
x" ( s )
C  sup
2
s( 0 , h )
12
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
x0  x(0)
xk 1  xk  h  (h, xk )
Il metodo di Eulero è un caso particolare:
(h, xk )  f (h, xk )
13
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 h0)
2. l’errore globale di approssimazione (cosa
succede se itero il procedimento N volte ed N∞)
14
1. Misura della qualità di un metodo
ad un passo
Errore locale di discretizzazione in x(0) :
e( x0 , h)  ( x0 , h)  ( x0 , h)
 x(h)  x(0)

, se h  0
( x0 , h)  
h
 f (0, x(0)), se h  0
15
Dimostro che
lim e( x(0), h)  0
h 0
Svolgimento:
x(h)  x(0)
lim e( x(0), h)  lim
 f (t , x(t )) 
h 0
h 0
h
x(h)  x(0)  hf (t , x(t ))
 lim

h 0
h
x(h)   x(0)  hf (t , x(t )) 
Ch
 lim
 lim
0
h 0
h 0
h
h
2
16
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.
17
2. Stima dell’errore globale di
discretizzazione
a. Illustrazione del problema tramite un esempio
b. stima
18
2.a Illustrazione del problema
tramite un esempio
Considero la seguente funzione:
Si verifica che soddisfa
una problema di Cauchy:
 x' (t )  x(t )

 x ( 0)  1
x(t )  e
t
x' (t )  e
x' (t )  x(t )
t
19
Sviluppo in serie di Taylor
 x' (t )  x(t )  f ( x(t ))

x(0)  1

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).
x1  x(h)  x(0)  h  x' (0)  1  h
Calcolo i punti successivi
x2  x1  h  f (1  h)  (1  h)  h  (1  h)  (1  h)


2
x3  x2  h  f 1  h  (1  h)  h  (1  h)  (1  h)
2
2
2
20
3
Quindi la funzione è:
e le iterate sono:
x(t )  e
t
x k  (1  h)
k
Voglio stimare l’errore (t=hk, così lo stimo nei punti della
successione tk)
kh
k
kh
k log(1 h )
k
x(kh)  x  e  (1  h)  e  e
21
Promemoria: Teorema di Lagrange
Sia f continua in [a,b], derivabile in (a,b). Allora esiste un punto c in (a,b) tale che
f (b)  f (a)  f ' (c)(b  a)
22
Usando il teorema di Lagrange
Con a=klog(1+h) e b=kh
k log(1 h )
  (log( 1  h), h)

e e
 e (kh  k log( 1  h)) 
 k (h  log( 1  h))  0
kh
Perché >0, quindi e >1
x(kh)  xk   
k  
Si può dimostrare anche che:
p
x(kh)  xk
 
p
k  
x
23
2.b Errore globale di
discretizzazione
Misura della qualità della approssimazione a
tempi molto maggiori di quello in cui è data la
condizione iniziale
E (h)  sup x(t )  xk
t[ 0 ,T ]
24
Stima di E(h) per il metodo di Eulero
Ipotesi ulteriori: f ed f ‘ limitate in R
M  sup f (t ) ; L  sup f ' (t ) ;
tR
tR
Fissato T e scelto h=T/N, N intero positivo fissato, vogliamo dimostrare che
Tesi
xN  x(T )  Mh(e
LT
1) / 2
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 h0
Osservazione 4: se T è grande, h dovrà essere scelto molto piccolo se si vuole avere
una buona approssimazione.
25
• Svolto fino a qui il 16 nov 09
26
M  sup f (t ) ; L  sup f ' (t ) ;
tR
tR
xN  x(T )  Mh(e LT 1) / 2
Dimostrazione
Osservazione: se x è soluzione del PC-EDO, allora
t
x(t )  x( s)   f ( , x( )) d  x( s) 
s
Caso particolare:
t s
 f (  s, x(  s))d ,
s, t
0
h
x(kh)  x(( k  1)h)   f (  (k  1)h, x(  (k  1)h)) d
0
Ponendo
d k  xk  x(kh)
si ha dunque
Ometto la
dipendenza da t
d k   yk 1  hf ( yk 1 ) 
h


  x(( k  1)h)   f (  (k  1)h, x(  (k  1)h)) d27
0


d k   yk 1  hf ( yk 1 ) 
h


  x(( k  1)h)   f (  (k  1)h, x(  (k  1)h)) d  
0


h
 d k 1   ( f ( xk 1 )  f ( x(( k  1)h   )) d
0
Applicando il teorema
di Lagrange e le ipotesi di limitatezza di f ed f ‘ si ha
h
d k  d k 1  L  ( x(( k  1)h   )  xk 1 d 
0
h
 d k 1  L  ( x(( k  1)h   )  x(( k  1)h)  x(( k  1)h)  xk 1 d 
0
h
 d k 1  L  ( x(( k  1)h   )  x(( k  1)h)  x(( k  1)h)  xk 1 d
0
Vd pagina successiva
|dk-1|
28
Osservazione: per t=(k-1)h+ e s=(k-1)h l’equazione
t
x(t )  x( s)   f ( , x( )) d  x( s) 
s
diventa
t s
 f (  s, x(  s))d ,
s, t
0

Ipotesi di limitatezza
x(( k  1)h   )  x(( k  1)h)   f ( x(( k  1)h   )) d  M
0
Pertanto
h
d k  d k 1  Lh d k 1  LM  d
0
d k  (1  Lh) d k 1
2
h
 LM , k
2
Iterando questa relazione ‘all’indietro’ fino a k=0 e
ricordando che d0=x(0)-y0=0 si trova:
29
Iterando questa relazione ‘all’indietro’ fino a k=0 e
ricordando che d0=x(0)-y0=0 si trova:
dk
 LMh 1  (1  Lh)  (1  Lh)
2
2
   (1  Lh)
k 1
/ 2 
 Mh[(1  Lh) k  1] / 2
Scegliendo ora k=n e ricordando che h=T/N,
osservando che (1+LT/N)N<eLT
xN  x(T )  Mh(e
LT
1) / 2
30
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
lim e(t , h)  0
h 0
e(t , h)
lim
c
p
h 0
h
Ovvero che l’errore locale
tende a zero con la stessa
rapidità di hp per h0
In tal caso si dice che l’errore
locale è un infinitesimo di ordine p
31
e(t,h)=O(hp)
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 è
x(t  h)  x(t )  hx' (t )  h2 x' ' (t ) / 2    h p x( p) (t ) / p! Rp (h)
Da cui si ricava
x(t  h)  x(t )
 x' (t )  hx' ' (t ) / 2    h p 1 x ( p ) (t ) / p! R p (h) / h
h
32
Poiché x è soluzione del PC-EDO
x' (t )  f ( x(t ))


x' ' (t )  f ' ( x(t )) x' (t )  f ' ( x(t )) f ( x(t ))

 x' ' ' (t )  f ' ' ( x(t )) f 2 ( x(t ))  ( f ' ( x(t ))) 2 f ( x(t ))

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
( x, h)  f ( x)  hf ' ( x) f ( x) / 2  h 2 [ f ' ' ( x) f 2 ( x)  ( f ' ( x)) 2 f ( x)] / 6
E la successione come
xk 1  xk  h( xk , h)
33
Calcolo l’errore locale di discretizzazione in x(0)
e( x0 , h)  ( x0 , h)  ( x0 , h)
:
 x(h)  x(0)

, se h  0
( x0 , h)  
h
 f (0, x(0)), se h  0
Ricordo anche l’espansione in serie di Taylor
x" (0)
3 x ' " ( 0)
x(h)  x(0)  h  x' (0)  h 
h 
 R3 (h)
2
3!
2
Da cui
R3 (h)
x" (0)
3 x ' " (0)
x(h)  x(0)  h  x' (0)  h 
h 

2
3!
h
2
quindi
R3 (h) 
x" (0)

2 x ' " (0)
e(t , h)   x' (0)  h
h 

  (t , h)

2
3!
h 

Sostituendo l’espressione di (xo,h) :
R3 (h)
e( x0 , h) 
h
3
34( h )
e( x0 , h)  O
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.
35
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.
36
Ipotesi:  derivabile due volte
( x0 , h)   ( x0 ,0)  h h ( x0 ,0)  O(h )
2
x" (0)
3 x ' " ( 0)
x(h)  x(0)  h  x' (0)  h 
h 
 R3 (h)
2
3!
2
f ' ( x0 ) f ( x0 )
e( x0, h)  f ( x0 )  h
 O( h 2 ) 
2
 ( x0 ,0)  h n ' ( x0 ,0)  O(h 2 )
Per ottenere un metodo
di ordine due basterà
allora scegliere in
modo tale che per ogni
x si abbia
=0
37
Esempio
( x0 , h)  Af ( x0 )  Bf ( x  Chf ( x0 ))
Con A, B, C parametri non negativi da determinare. La condizione
f ' ( x0 ) f ( x0 )
f ( x0 )  h
 ( x0 ,0)  h n ' ( x0 ,0)  0
2
diventa
f ' ( x0 )
f ( x0 )  h
 ( A  B) f ( x0 )  hBCf ' ( x0 ) f ( x)
2
Che è verificata se A, B, C soddisfano il sistema di
equazioni
A

B

1


 BC  1

2

38
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
39
Metodo di Runge-Kutta
Si ottiene definendo
( x0, h)  (k1  2k2  2k3  k4 ) / 6
Con
Si può dimostrare che è un metodo ad
un passo di ordine 4.
k1  f ( x0 )
 k  f ( x  hk / 2)
 2
0
1

 k3  f ( x0  hk 2 / 2)
 k 4  f ( x0  hk3 )
40
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 modificato
E(h)=O(h4) per il metodo di Runge-Kutta
41
Nota sui sistemi del primo ordine
I metodi visti finora si estendono ad un
sistema di N equazioni del primo ordine.
Che in forma vettoriale diventa
x' (t )  f ( x(t ))
x(0)  x0
Lo spazio euclideo RN dove è assegnata la
condizione iniziale e dove evolve la traiettoria x(t)
prende il nome di spazio delle fasi
 x'1 (t )  f1 ( x(t ))
 x' (t )  f ( x(t ))
2
 2



 x' N (t )  f N ( x(t ))

x1 (0)  x1

 x2 (0)  x2



 x (0)  x
N
N

42