LUISS
Corso di
Metodi matematici per economia e finanza.
Prof. F. Gozzi
a.a. 2009/10
Integrazione numerica di equazioni differenziali
- Prof.ssa G. Rotundo
1
Elenco degli argomenti
Formule di integrazione numerica di EDO 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 Heun
3. Metodo di Runge-Kutta
Formule per il calcolo di integrali
Metodo dei rettangoli
Metodo dei trapezi
Metodo di Simpson 1/3
Corrispondenza tra metodi di calcolo numerico di EDO ed integrali
Corrispondenza tra il metodo di Eulero (modificato) ed il metodo dei rettangoli
Corrispondenza tra il metodo di Heun ed il metodo dei trapezi
Modello preda-predatore
Descrizione al modello
Formalizzazione tramite un sistema di equazioni differenziali
Linearizzazione
Integrale primo, grafici di curve di livello e traiettorie.
2
Testi di riferimento
MATLAB – manuale di riferimento
R. A. Adams, Calcolo Differenziale 2, Casa Editrice
Ambrosiana, ISBN 978 – 88 – 408 – 1390 -5
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
16 nov 09
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 4
alcuni metodi numerici per approssimare la soluzione.
Metodi numerici di integrazione di EDO
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 Heun
2. Metodo di Eulero modificato
3. 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:
x(t  h)  x(t )
lim
 x' (t )
h 0
h
Per h abbastanza piccolo il rapporto incrementale è quasi uguale
alla derivata:
x(t  h)  x(t )
 x' (t )
Sostituendo
si ha che
da cui
h
x' (t )  f (t , x(t ))
x(t  h)  x(t )
 f (t , x(t ))
h
x(t  h)  x(t )  h  f (t , x(t ))
7
La relazione
x(t  h)  x(t )  h  f (t , x(t ))
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

8
x(t  h)  x(t )  h  f (t , x(t ))
In corrispondenza a tk si può definire una successione di xk :
t0
t1
t2

t k 1
x0  x(t0 )
x1  x(t1 )  x(t0  h)  x0  h  f (t0 , x0 )
x2  x(t 2 )  x(t1  h)  x1  h  f (t1 , x1 )

xk 1  x(t k 1 )  x(t k  h)  xk  h  f (t k , xk )
9
La funzione approssimante, lineare a tratti, può essere descritta, t: 0tT, dalla
espressione:
xh (t )  xk  (t  hk ) f (t , xk ), t [kh, (k  1)h), k  0,1,2...
x2
x1
x0
t0
t1
t2
Interpretazione geometrica del metodo di Eulero: il punto successivo è ottenuto
dal precedente spostandosi con pendenza fissata per un intervallo di ampiezza
h. La pendenza è data dalla funzione f calcolata nel primo estremo di ciascun
10
intervallo.
Esempio 1
Considero la seguente funzione:
x(t )  t
2
Si verifica che soddisfa una equazione differenziale, infatti:
x' (t )  2t
t  x' (t )  2t 2
t  x' (t )  2 x(t )
x' (t )  2 x(t ) / t
Quindi
x(t )  t
2
è soluzione della equazione differenziale:
x' (t )  f (t , x(t ))
f (t , x(t ))  2 x(t ) / t
con
11
Grafico della soluzione x(t)=t2 e della sua
approssimazione tramite il metodo di Eulero
Considero la condizione iniziale
t0=0.5 e x(0.5)=0.25;
passo h=0.2: la generica iterazione
col metodo di Eulero è:
xk 1  xk  0.2  2 xk / tk 
6
5
4
x(t  h)  x(t )  h  f (t , x(t ))
3
2
Problema: la successione
approssimante (spezzata in nero con
i punti evidenziati) si allontana ben
presto dalla soluzione .
1
0
0.5
1
1.5
2
2.5
E’ necessario dare una stima della bontà della approssimazione:
x(t  h)  x(t )  h  f (t , x(t ))
Ricordando che:
x’(t)=f(t,x(t))
12
Confronto con lo 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
da cui ricavo una stima dell’errore dell’approssimazione di x(h) tramite x(0)+hx’(0) :
x" ( s)
2
x(h)  x(0)  h  x' (0)   h 
 Ch
2
2
Ammesso che esista
x" ( s )
C  sup
2
s( 0 , h )
13
(*)
Esempio 2
Considero la seguente funzione:
Si verifica che soddisfa
un problema di Cauchy:
 x' (t )  x(t )

 x ( 0)  1
x(t )  e
t
x' (t )  e
x' (t )  x(t )
t
14
Costruzione della approssimazione
numerica della soluzione di
Osservo che la variabile t non appare
da sola, ma soltanto tramite x, quindi
l’approssimazione di Eulero diventa:
 x' (t )  x(t )  f ( x(t ))

x(0)  1

xk 1  xk  h  f ( xk )
Osservo che la funzione f(.) è la funzione identica (restituisce l’argomento, immutato).
Quindi l’approssimazione di Eulero diventa:
xk 1  xk  h  xk
Scrivo il primo punto che si ottiene tramite l’approssimazione con il metodo di Eulero:
x1  x(h)  x(0)  h  x' (0)  1  h  x(0)  1  h
Scrivo i punti successivi che si ottengono applicando il metodo:
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
15
3
x(t )  e
Quindi la funzione soluzione del problema è:
t
x k  (1  h)
e le iterate approssimanti sono:
5
x(t )  e
4.5
k
t
t
x(t)
xk
0 0.00 1.00
1.00
1 0.40 1.49
1.40
3.5
2 0.80 2.23
1.96
3
3 1.20 3.32
2.74
2.5
4 1.60 4.95
3.84
x(t)
x4 4
x3
x2
2
x11.5
x0 10
0.2
0.4
0.6
0.8
t
1
1.2
1.4
1.6
Osservo che anche in
questo caso la distanza
(barre verticali in nero)
della soluzione da quella
16
approssimata aumenta al
crescere di k.
k
Programma matlab eulero.m
h=0.4;
T=h*4;
t=[0:h:T]; %t assume 5 valori: 0, 0.3, 0.6, 0.9, 0.12
%gli indici slittano di una unità per via della memorizzazione in MATLAB
for i=1:1:length(t)
xk(i)=(1+h)^(i-1);
end
x=exp(t);
plot(t,x,'k',t,xk,'r+-');xlabel('t');ylabel('x(t)');
for i=1:1:length(t)
fprintf(' %2d | %4.2f | %4.2f | %4.2f \n',i,i, t(i),x(i),xk(i));
end
17
Nota
Da questo momento in poi consideriamo, per semplicità, il problema
autonomo
 x' (t )  f ( x(t ))

 x(t0 )  x0  X
t  T
18
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  ( xk , h)
Il metodo di Eulero corrisponde ad una particolare scelta di :
( xk , h)  f ( xk )
19
Analisi della qualità del metodo
L’errore di approssimazione è piccolo se h è piccolo?
Voglio determinare
1. se l’errore di approssimazione è piccolo
se h è piccolo
(procedimento:
(a) definisco una misura per l’errore e
(b) dimostro che il limite è 0 per h0)
2. l’errore globale di approssimazione
(devo esaminare cosa succede se itero
l’approssimazione N volte ed N∞)
20
1. (a)Misura della qualità di un
metodo ad un passo
Definisco l’errore locale di discretizzazione in x(0) :
e( x0 , h)  ( x0 , h)  ( x0 , h)
con
 x(h)  x(0)

, se h  0
( x0 , h)  
h
 f ( x(0)), se h  0
21
e( x(0), h)  0
1.(b) Dimostro che lim
h 0
Svolgimento:
x ( h )  x ( 0)
lim e( x(0), h)  lim
 f ( x(t )) 
h 0
h 0
h
x(h)  x(0)  hf ( x(t ))
 lim

h 0
h
utilizzo (*)
x(h)   x(0)  hf ( x(t )) 
Ch
 lim
 lim
0
h 0
h 0
h
h
2
22
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
nel paragrafo che segue.
23
2. Stima dell’errore globale di
discretizzazione
(a)Illustrazione del problema tramite un esempio:
dipendenza della qualità della approssimazione dal
numero di iterazioni.
(b) Stima dell’errore globale di discretizzazione:
dipendenza dall’ampiezza del passo.
24
2.(a) Illustrazione del problema tramite un esempio
 x' (t )  x(t )

 x ( 0)  1
Considero nuovamente il problema di Cauchy
La funzione soluzione del problema è:
x(t )  e
e le iterate della approssimazione di Eulero sono:
x k  (1  h)
t
Voglio stimare l’errore, cioè la differenza tra la soluzione esatta e quella approssimata.
Pongo t=hk, così stimo l’errore solo nei punti della successione tk .
x(tk )  xk  x(kh)  xk  e  (1  h) 
kh
k log(1 h )

e

e
Devo semplificare questa espressione in modo da capirne facilmente
kh
k
l’andamento quando il numero k di passi aumenta (all’infinito).
25
k
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)
La dimostrazione fa uso del teorema di Rolle, che a sua volta fa uso del
teorema di Weierstrass e sono riportati sul libro di matematica generale.
26
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
27
Tramite l’esempio ho dimostrato che la qualità
dell’approssimazione peggiora (la distanza tra
soluzione esatta e soluzione approssimata va
all’infinito) se il numero delle iterazioni aumenta.
Un altro aspetto del problema riguarda la dipendenza
della qualità della approssimazione dal passo h.
Nel prossimo paragrafo dimostramo, in generale, che
se il passo h diminuisce la qualità della
approssimazione migliora.
28
2.b Errore globale di
discretizzazione
Misura della qualità della approssimazione a
tempi molto maggiori di quello in cui è data la
condizione iniziale
Definizione: l’errore globale della approssimazione è dato da:
E (h)  sup x(t )  xk
t[ 0 ,T ]
29
23 nov 09
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.
30
xN  x(T )  Mh(e LT 1) / 2
M  sup f (t ) ; L  sup f ' (t ) ;
tR
tR
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 ( x(  s))d ,
s, t
0
h
x(kh)  x(( k  1)h)   f ( x(  (k  1)h)) d
Ponendo
d k  xk  x(kh)
0
si ha dunque
d k  xk 1  hf ( xk 1 ) 
h


  x(( k  1)h)   f ( x(  (k  1)h)) d 
0


31
d k  xk 1  hf ( xk 1 ) 
h


  x(( k  1)h)   f ( 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
ora sommo e sottraggo la stessa quantità
 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|
32
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 ( 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:
33
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
q.e.d
34
Metodi di ordine superiore basati
sullo sviluppo di Taylor
Vogliamo ora studiare dei metodi di approssimazione che, a
parità di passo 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
35
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
36
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)
37
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 ( 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
x" (0)
3 x ' " (0)
x(h)  x(0)  h  x' (0)  h 
h 
 R3 (h)
2
3!
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 (x0,h) :
R3 (h)
e( x0 , h) 
h
3
38( h )
e( x0 , h)  O
Osservazioni
1. Questi metodi possono essere usati per definire
metodi di approssimazione ad un passo di ordine
comunque 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.
39
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.
40
Ipotesi:  derivabile due volte
( x0 , h)   ( x0 ,0)  h 'h ( x0 ,0)  O(h )
2
x" (0)
x(h)  x(0)  h  x' (0)  h 
 R2 (h)
2
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
41
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 )
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

42
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
Corrisponde ad utilizzare come pendenza la media delle pendenze
all’inizio ed alla fine di ciascun intervallo della partizione su t.
Metodo di Eulero modificato: A=0, B=1, C=1/2
Corrisponde ad utilizzare come pendenza la pendenza calcolata nel punto
intermedio di ciascun intervallo della partizione.
43
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 )
44
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 può
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
45
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

46
30 nov 09
Metodi numerici per il calcolo
approssimato di integrali
•Metodo dei rettangoli
•Metodo dei trapezi
•Metodo di Simpson
47
Osservazione
Nel caso in cui f(t,x(t)) dipende solo da t il problema
di Cauchy corrispondente si riduce al calcolo di
una particolare primitiva di f(t).
 x' (t )  f (t ) t  T

 x(t0 )  x0  X
t
x(t )   f ( )d  x0
t0
48
Obbiettivo: calcolare l’integrale di
una funzione reale su [a,b]
Strumenti noti da matematica generale: il teorema
fondamentale del calcolo integrale ( teorema di
Newton-Leibniz), che risolve il problema nel caso
in cui f è continua ed è nota una primitiva di f,
cioè una funzione derivabile tale che F’(t)=f(t)
per ogni x: a<x<b. In tal caso vale la formula:
b
 f (t )dt  F (b)  F (a)
a
49
Perché cambiare metodo?
L’applicabilità pratica di questo metodo analitico è
limitata, infatti:
1. Esistono funzioni continue che si incontrano
frequentemente nelle applicazioni dell’analisi
matematica la cui primitiva non può essere
espressa in termini di funzioni elementari ( per
esempio sen(t)/t, 1/((1-t2)(1-k2t2))1/2, exp(-t2);
2. Il calcolo può essere complicato ( p.es. nel
caso di integrazione di funzioni razionali fratte,
che si basano sulla conoscenza di radici di
polinomi).
50
Obbiettivo: calcolare
UN VALORE APPROSSIMATO
dell’integrale
b
 f (t )dt
a
Strumenti: metodi di integrazione numerica facilmente
implementabili:
Rettangoli, Trapezi, Simpson [...]
Idea: integrare, invece di f, una sua approssimazione g che si
sappia integrare in maniera esatta (a meno, ovviamente, di
errori di arrotondamento).
51
Punto fondamentale
Calcolare a priori l’errore commesso con la specifica
procedura adottata. Serve per garantire l’affidabilità del
metodo e bisognerà pertanto stimare
E
b
b
a
a
 f (t )dt   g (t )dt
In ciascun metodo la costruzione di g è fatta definendo una partizione
dell’intervallo [a,b]. Precisamente, si fissa un intero positivo N e si divide
l’intervallo di integrazione [a,b[ in parti uguali con i punti
tk  a  k (b  a) / N , k  0,1,, N
Il numero h=(b-a)/N è il passo della discretizzazione
52
METODO DEI RETTANGOLI
a
b
Devo calcolare l’area sottesa dal grafico della funzione.
53
METODO DEI RETTANGOLI
Idea: costruisco la partizione di [a,b]
tk  a  k (b  a) / N , k  0,1,, N
Obbiettivo: approssimare l’area con rettangoli di base (tk+1-tk) ed altezza uguale
all’altezza della funzione nel punto intermedio di ciascun intervallo (tk+tk+1)/2
Preparo gli elementi che
servono per questa
approssimazione
a=t0
t1
t2
b=tN
54
Considero i punti intermedi in ciascun intervallo
a=t0
t1
t2
b=tN
55
Considero il valore della funzione in ciascun punto intermedio in ciascun
intervallo
a=t0
t1
t2
b=tN
56
Considero la funzione costante a tratti, avente, su ciascun intervallo, il
valore costante uguale all’altezza della funzione calcolata nel punto
intermedio dell’intervallo
 t k  t k 1 
g (t )  f 
, t  [t k , t k 1 ), k  1,2, , N
 2 
a=t0
t1
t2
b=tN
57
Considero la funzione costante a tratti, avente, su ciascun intervallo, il
valore costante uguale all’altezza della funzione calcolata nel punto
intermedio dell’intervallo
 t k  t k 1 
g (t )  f 
, t  [t k , t k 1 ), k  1,2, , N
 2 
Ovviamente g coincide con f in tutti i punti intermedi di ciascun intervallo.
a=t0
t1
t2
b=tN
58
Approssimo l’area da calcolare con la somma delle aree (base x altezza) di questi
rettangoli. Questo calcolo è immediato:
•tutte le basi hanno ampiezza costante (b-a)/N
•le altezze sono date dal valore della funzione nel punto intermedio dell’intervallo
N 1
I 
k 0
 t k  t k 1  b  a
f

 2  N
a=t0
t1
t3
b=tN
59
Osservazioni
• Il metodo dei rettangoli è ispirato dalla
definizione di integrale definito.
• I è una particolare somma integrale.
• Il numero I dipende dalla partizione scelta
e quindi dal passo h (h=(b-a)/N).
60
Stima dell’errore
E
b
b
a
a
 f (t )dt   g (t )dt
61
MEMO: da matematica generale
Proprietà dell’integrale definito
L’area si può calcolare dividendola in due parti: (a, c) e (c, b)
b
c
b

a
f ( x)dx   f ( x )dx   f ( x )dx
ya
0
c
f(x)
x
a
c
b
62
Applico la proprietà dell’integrale ripetutamente:
b

a
t1
t2
a
t1
f (t )dt   f (t )dt   f (t )dt   
b

f (t )dt 
t N 1
N 1 t k 1


f (t )dt
k 0 tk
Posso ripetere il procedimento per g e quindi
E
b
b
N 1 t k 1
a
a
k 0 tk
 f (t )dt   g (t )dt     f (t )  g (t )dt
63
MEMO
ab  a  b
In generale
a  b  c  d   a  b  c  d 
Il valore assoluto di una somma è minore od uguale alla somma dei valori assoluti
64
b
b
N 1 t k 1
a
a
k 0 tk
 f (t )dt   g (t )dt     f (t )  g (t )dt 
E
N 1 t k 1


k 0 tk
N 1 t k 1
f (t )  g (t ) dt  

k 0 t k
 t k  t k 1 
f (t )  f 
 dt 
 2 
Scrivo l’espressione di g
Ho applicato la proprietà del valore assoluto
MEMO:
teorema di Lagrange: sotto opportune ipotesi f(b)-f(a)=f ‘ (c ) (b-a), a<c<b
Applico il teorema considerando b=t e a= (tk+tk+1)/2, c=k , t< k < (tk+tk+1)/2
N 1 t k 1


k 0 t k
t k  t k 1
f ' ( k ) t 
dt
2
65
Osservazione
t
tk
(tk+tk+1)/2
tk+1
t k  t k 1 h
t
 , t  [t k , t k 1 )
2
2
Un punto t dista dal centro dell’intervallo meno di metà della lunghezza
dell’intervallo.
66
Utilizzo questa proprietà
per svolgere il passaggio
N 1 t k 1


k 0 tk
N 1 t k 1


k 0 tk
t k  t k 1 h
t
 , t  [t k , t k 1 )
2
2
t k  t k 1
f ' ( k ) t 
dt  
2
k 0
N 1 t k 1

tk
h
f ' ( k ) dt 
2
N 1 t k 1
h
h
sup f ' ( k ) dt  sup f ' ( k )   1dt 
2
2 k 0 t k
a  b
a  b
h  N 1
h

 sup f ' ( k )   (t k 1  t k )   sup f ' ( k ) ( Nh) 
2  k 0
2
a  b
 a  b
h
h
 sup f ' ( k ) ( Nh)  sup f ' ( k ) (b  a)
2
2
a  b
a  b
67
Conclusione: ho dimostrato che
h
E (h)  sup f ' ( k ) (b  a)
2
a  b
68
METODO DEI TRAPEZI (o di Crank-Nicolson)
Mentre nel metodo dei rettangoli
si approssima f con g costante
a tratti,
nel metodo dei trapezi si
approssima f con g
costruita partendo dalle
rette secanti in ciascun
intervallo
a
=
t0
t1
a
=
t0
t1
t2
b=t
N
Quindi la funzione g è lineare in ogni intervallo
e coincide con f negli estremi di ciascun
intervallo.
Il metodo si chiama ‘dei trapezi’
perché approssimo l’area tramite
trapezi (eventualmente degeneri)
invece che tramite rettangoli.
t2
b=t
N
69
L’approssimazione corrisponde al calcolo dell’area evidenziata in
rosso.
Osservazione:
alcuni trapezi possono essere degeneri, per esempio
composti da triangoli o rettangoli
Obbiettivi:
•calcolo dell’area
•stima dell’errore
70
Osservazione: rispetto alla normale visualizzazione, i
trapezi sono ruotati di 90 gradi:
Base minore
altezza
Base maggiore
MEMO: area di un trapezio
(base minore + base maggiore)*altezza /2
71
MEMO: da matematica generale
Equazione di una retta che passa per i punti
y
y2
(x1,y1) e (x2,y2)
y2  y1
m
x2  x1
y2  y1
y  y1
m
x2  x1
x  x1
y  y1 y2  y1

x  x1 x2  x1
y2  y1
y  y1 
( x  x1 )
x2  x1
y
y1
da cui
x
da cui
da cui
x1
x
x2
y2  y1
y
( x  x1 )  y1
x2  x1
72
Calcolo della retta secante
nel singolo intervallo
Applico la formula
y2  y1
y
( x  x1 )  y1
x2  x1
Considerando le particolari
coordinate dei punti:
f (t k 1 )  f (t k )
g (t ) 
(t  t k )  f (t k )
t k 1  t k
f(tk+1)
f(t)
f(tk)
…
tk
t
tk+1
…
t  [t k , t k 1 ), k  0,1,..., N  1
Osservo che la funzione g coincide con la funzione f negli estremi dei
singoli intervalli ed è lineare in ogni intervallo.
73
t
Area approssimata = somma delle
aree dei singoli trapezi
Base minore+base maggiore
b

altezza
N 1
f (t )dt    f (t k )  f (t k 1 )   (t k 1  t k ) / 2 
k 0
a
ba
   f (t k )  f (t k 1 )   
/2
 N 
k 0
N 1
Perché gli intervalli hanno tutti la stessa ampiezza
74
Stima dell’errore
Ipotesi in più: f derivabile due volte in (a,b) con derivata f’’ limitata.
sup f ' ' ( )  M
 ( a ,b )
Calcolo l’errore nel singolo intervallo, poi sommo per ottenere l’errore
totale.
Inizio il calcolo partendo dall’espressione di g
f (t k 1 )  f (t k )
g (t ) 
(t  t k )  f (t k )
t k 1  t k
t  [t k , t k 1 ), k  0,1,..., N  1
75
Per il teorema di Lagrange
f (t )  f (tk )  f ' ( k )(t  tk )
Inoltre g ed f coincidono sui punti della partizione
f (t k 1 )  f (t k )
g (t )  g (t k )  g (t )  f (t k ) 
(t  t k )
t k 1  t k
Sottraggo

f (tk 1 )  f (tk ) 
(t  tk )
( f (t )  f (tk ))  ( g (t )  g (tk ))   f ' ( k ) 
tk 1  tk


Da cui

f (t k 1 )  f (t k ) 
 (t  t k )
f (t )  g (t )  ( f (t )  f (t k ))  ( g (t )  g (t k ))   f ' ( k ) 
t k 1  t k


Ora cerco di semplificare questa quantità
76

f (t k 1 )  f (t k ) 
 f ' ( k ) 

t k 1  t k


Considero quindi
Lavoro dapprima su
f (tk 1 )  f (tk )
f (tk 1 )  f (tk )  f ' (k )tk 1  tk 
Applico la formula di Lagrange ed ottengo
Sostituisco ed ottengo:

f (t k 1 )  f (t k ) 
 f ' ( k ) 
  f ' ( k )  f ' ( k ) t k 1  t k
t k 1  t k


77

f (t k 1 )  f (t k ) 
 f ' ( k ) 
  f ' ( k )  f ' ( k ) t k 1  t k
t k 1  t k


Applico il teorema di Lagrange alla funzione derivata prima:
f ' ( k )  f ' ( k )  f ' ' ( )( k   k ) 
 f ' ' ( )t k 1  t k   f ' ' ( ) t k 1  t k
Sostituisco ed ottengo che in ciascun intervallo

f (tk 1 )  f (tk ) 
 tk 1  t k 
f (t )  g (t )   f ' ( k ) 
tk 1  t k


 f ' ' ( ) t k 1  tk   f ' ' ( ) h  sup f ' ' ( ) h
2
2
 ( a ,b )
2
78
Errore totale: sommo sugli N
intervalli
b
E
b
b
 f (t )dt   g (t )dt  
a
N 1 t k 1


k 0 tk
a
f (t )  g (t ) dt 
a
N 1 t k 1
f (t )  g (t ) dt   Mh dt  Mh
2
k 0 tk
N 1
N 1
k 0
k 0
t k 1
N 1
k 0
2
 dt 
tk
  Mh 2 (t k 1  t k )  Mh 3 NMh 3
Quindi la stima cercata è:
E  E (h)  NMh3  Mh 2 (b  a)
79
METODO DI SIMPSON
I passi effettuati finora per il calcolo approssimato
dell’integrale hanno portato da una prima
approssimazione mediante una funzione
costante a tratti ad una approssimazione tramite
una funzione rettilinea a tratti, in cui l’errore va a
zero più velocemente quando h0.
Il passo successivo riguarda l’approssimazione
mediante una funzione quadratica.
Passi ulteriori possibili riguardano
l’approssimazione della funzione assegnata
mediante polinomi di ordine maggiore.
80
Interpolata quadratrica
Considero una funzione in [-1,1] e la sua interpolata
quadratica definita da:
g (t )  At  Bt  C , t [1,1]
2
Dove le tre costanti sono determinate dalle tre condizioni
g (1)  A  B  C  f (1)
g (0)  C  f (0)
g (1)  A  B  C  f (1)
che chiedono che le
funzioni coincidano su
quei tre punti.
Osservo che l’unica soluzione del sistema lineare è:
f (1)  f (1)
f (1)  f (1)
A
 f (0), B 
, C  f (0)
2
2
81
Interpolata quadratrica
 t

t
 g (t )dt   At  Bt  C dt  A 3  B 2  Ct 
1
Quindi
1
3
2
1
2
1
1

 1
A B
A B
2
   C  [   C ]  A  2C 
3 2
3 2
3

f (1)  4 f (0)  f (1)

3
Questo calcolo si generalizza facilmente ad un arbitrario intervallo [c,d] tramite un
cambio di variabile
cd d c
t

t
2
2
Che
trasforma l’intervallo [-1,1] in [c,d]. L’integrale diventa quindi
d
2
c g (t )dt  3 A  2C  d  c f (c)  4 f ((c  d ) / 2)  f (d )/ 6
82

Su un intervallo qualsiasi [c,d]
d
h


g
(
t
)
dt

f
(
c
)

4
f
(
m
)

f
(
d
)
c
3
g(x)
f(x)
con semiampiezza dell’intervallo h
h=(m-c)=(d-m)
c
h
m
h
d
Il nome con cui il metodo è più conosciuto è “Simpson 1/3”, per distinguerlo da un altro
metodo di Simpson in cui la costante è diversa (“Simpson 3/8”)
83
Su più intervalli consecutivi
Applico la formula considerando i punti della partizione a gruppi di tre:
iniziale ( c ), finale (d) ed intermedio (m) e sommando tutti i risultati
t0
t1
t2
Quindi la funzione g è lineare in ogni intervallo e coincide con f negli
estremi di ciascun intervallo.
b=tN
84
Su più intervalli consecutivi
Applico la formula considerando i punti della partizione a gruppi di tre: iniziale ( c ),
finale (d) ed intermedio (m) e sommando tutti i risultati.
b

a
1  b  a  N 2
f (t )dt  
   f (t k )  4 f (tk 1 )  f (t k  2 ) 
3  N k 0
k pari
(h=(b-a)/N)
h
t0
h
t1
h
t2
h
h
h
b=tN
85
Osservazione
Le formule ottenute per la quadratura sono
casi particolari di formule di Newton-Cotes
in quanto la partizione sull'asse delle
ascisse è costituita da punti equidistanti.
Le formule si dicono di tipo chiuso quando
x0=a, xn=b (per esempio metodo dei trapezi
e metodo di Simpson);
si dicono di tipo aperto quando x0=a+q, b=b-q
(per esempio nel metodo dei rettangoli).
86
Osservazione
Nel caso in cui f(t,x(t)) dipende solo da t il problema di Cauchy
corrispondente si riduce al calcolo di una particolare
primitiva di f(t) .
In questo caso si può dimostrare che il metodo di Heun si
riduce al metodo dei trapezi per il calcolo di

 f (t )dt
Si può anche dimostrare0che il metodo di Eulero modificato si
riduce a quello dei rettangoli, il metodo di Runge-Kutta a
quello di Simpson.
Nella prossima sezione mostriamo la corrispondenza tra
– il metodo di Eulero (modificato) ed il metodo dei rettangoli;
– il metodo di Heun ed il metodo dei trapezi.
87
14 dic 09
Confronto tra metodi per la risoluzione
di equazioni differenziali e metodi per il
calcolo di integrali
Mostriamo la corrispondenza tra
• il metodo di Eulero (modificato) ed il metodo dei
rettangoli;
• il metodo di Heun ed il metodo dei trapezi.
88
Confronto tra metodi per la risoluzione di EDO
e metodi per il calcolo di integrali
Formule di integrazione numerica di EDO
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 Heun
3. Metodo di Runge-Kutta
Formule per il calcolo di integrali
Metodo dei rettangoli
Metodo dei trapezi
Metodo di Simpson 1/3
 x' (t )  f (t , x(t ))

 x(t0 )  x0  X
t  T
Nel caso in cui
f(t,x(t)) dipende
solo da t il
problema di
Cauchy
corrispondente si
riduce al calcolo
di una particolare
primitiva di f(t).
 x' (t )  f (t ) t  T

 x(t0 )  x0  X
Mostriamo la corrispondenza tra il metodo di Eulero (modificato) ed il metodo dei
rettangoli e tra il metodo di Heun ed il metodo dei trapezi.
89
Corrispondenza tra il metodo di Eulero (modificato) ed il
metodo dei rettangoli
Procedimento: scrivo le iterazioni di Eulero, raggruppo i termini in maniera diversa,
mostro che, modificata, la formula è la stessa dei rettangoli.
Ricordo che il metodo di Eulero individua la successione approssimante è
definita per ricorrenza da
t0
t1
t2

t k 1
x0  x(t0 )
x1  x(t1 )  x(t0  h)  x0  h  f (t0 , x0 )
x2  x(t 2 )  x(t1  h)  x1  h  f (t1 , x1 )

xk 1  x(t k 1 )  x(t k  h)  xk  h  f (t k , xk )
90
scrivo xk+1 raggruppando i termini in modo diverso
xk 1  x(tk 1 )  x(tk  h)  xk  h  f (tk , xk )
x0  x(t0 )
x1  x0  h  f (t0 , x0 )
x2  x1  h  f (t1 , x1 )  ( x0  h  f (t0 , x0 ))  h  f (t1 , x1 ) 
 x0  h  f (t0 , x0 )  h  f (t1 , x1 )  x0  h   f (t0 , x0 )  f (t1 , x1 )  
1
 x0  h   f (tl , xl )
l 0
x3  x2  h  f (t 2 , x2 )  ( x1  h  f (t1 , x1 ))  h  f (t 2 , x2 ) 
 x0  h  f (t0 , x0 )  h  f (t1 , x1 )  h  f (t 2 , x2 ) 
2
 x0  h   f (t0 , x0 )  f (t1 , x1 )  f (t 2 , x2 )   x0  h   f (tl , xl )
l 0
(**)
quindi
k
xk 1  xk  h  f (t k , xk )  x0  h   f (tl , xl )
l 0
91
Quindi la trasformazione della approssimazione che si ottiene con
il metodo di Eulero
xk 1  x(tk 1 )  x(tk  h)  xk  h  f (tk , xk )
N 1
permette di scrivere
x N  x0  h   f (tl , xl )
l 0
Considero, per semplicità, il caso in cui x0=0 : x N  h 
N 1
 f (t , x )
l 0
l
l
N 1
Ricordo l’approssimazione dell’integrale che
 t k  t k 1 
I h f

si ottiene dal metodo dei rettangoli:

k 0

2
 t k  t k 1 
Le due formule coincidono se tl  
 , e questo di fatto
 2 
significa che stiamo usando il metodo di Eulero modificato.
92

Infatti, la modifica che viene considerata nel metodo di Eulero modificato
corrisponde a scrivere il generico passo della iterazione
xk 1  xk  h( xk , h)
con
( xk , h)  Af ( xk )  Bf ( xk  Chf ( xk ))
h
ponendo A=0, B=1, C=1/2, quindi  ( xk , h)  f ( xk  f ( xk ))
2
Reintroducendo esplicitamente la dipendenza dal tempo si ha:
h
h
(t k , xk , h)  f (t k  , xk  f (t k , xk ))
2
2
Poiché ora stiamo trattando il caso di funzioni che dipendono soltanto dal
tempo l’espressione si semplifica e si ha:  (t , x , h)  f (t  h )
k
k
k
2
Questo significa che xk+1 puo’ essere scritto, analogamente a quanto svolto in
k
(**), come
xk 1  xk  h   (t k , xk , h)  x0  h    (tl , xl , h)
l 0
Avendo posto, per semplicità, x0=0 , fissando k=N e sostituendo  si ottiene:
N 1
xN  h  
h
f (tl  )
2
Ovvero la stessa approssimazione del metodo dei rettangoli,
l 0
che considera come altezza del rettangolo il valore della
funzione nel punto intermedio di ciascun intervallo della partizione. Senza questa
modifica, la formula di Eulero considera rettangoli con altezza data dal valore della93
funzione nel primo estremo di ciascun intervallo.
Corrispondenza tra il metodo di Heun ed il metodo dei
trapezi
Nel metodo di Heun l’iterazione è definita da: xk 1  xk  h( xk , h)
con
( xk , h)  Af ( xk )  Bf ( xk  Chf ( xk )) con A=1/2, B=1/2, C=1
e quindi diventa
1
1
 ( xk , h)  f ( xk )  f ( xk  1 hf ( xk ))
2
2
Reintroducendo esplicitamente la dipendenza dal tempo si ha:
1
1
(t k , xk , h)  f (t k , xk )  f (t k  h, xk  hf ( xk ))
2
2
Poiché ora stiamo trattando il caso di funzioni che dipendono soltanto dal
tempo l’espressione si semplifica e si ha:
 (t k , xk , h) 
f (t k )  f (t k  h)
1
1
f (t k )  f (t k  h) 
2
2
2
ovvero il valor medio di f rispetto agli estremi del singolo intervallo della
partizione.
94
Avendo posto, per semplicità, x0=0 , fissando k=N e sostituendo  si ottiene:
N 1
xN  h  
l 0
f (tl )  f (tl  h)
2
Considerare come pendenza la media delle pendenze
agli estremi di ciascun intervallo della partizione
dovrebbe migliorare l’approssimazione rispetto al
metodo di Eulero, in cui si considera la tangente in uno
soltanto degli estremi del segmento.
t
t+h
Ricordo il valore che si ottiene con il metodo dei trapezi:
N 1
b
Base minore+base
maggiore
altezza
 f (t )dt    f (t
k 0
a
k
)  f (t k 1 )   (t k 1  t k ) / 2 
N 1
 f (t k )  f (t k 1 ) 
ba
   f (t k )  f (t k 1 )   
h
 / 2  
2
 N 

k 0
k 0 
N 1
Cioè esattamente ciò che si ottiene con il metodo di Heun
95
21 dic 09
Esempi numerici (in excel)
Mediante il metodo di
Eulero determinare i
valori approssimati
della soluzione del
problema ai valori
iniziali, nell’intervallo
[0,1], usando
(a) 5 passi di lunghezza
h=0.2
(b) 10 passi di
lunghezza 0.1
 dy
  x y
 dx
 y (0)  1
Calcolare l’errore ad ogni
passo, dal momento che il
problema coinvolge
un’equazione (lineare)
risolvibile in modo esatto la cui
soluzione è y=x-1+exp(-x)
96
Soluzione
(a) In questo caso abbiamo f(x,y)=x-y, x0=0,y0=1 e h=0.2, per cui xn=n/5,
yn+1=yn+0,2*(xn-yn) e l’errore è y-yn, per n=0,1,2,3,4,5. I risultati del
calcolo si ottengono tramite il foglio elettronico.
(b) h=0.1, n=10.
File excel: esercizio.xls
Osservazione
Si noti che l’errore alla fine del primo passo è circa un quarto dell’errore alla
fine del primo passo trovato nella partte (a), ma l’errore finale in x=1 è
solo la metà di quello della parte (a). Questo comportamento è
caratteristico del metodo di Eulero.
Se diminuiamo la lunghezza del passo h occorrono più passi per
raggiungere dal punto di partenza il valore particolare cui vogliamo
arrivare.
97
Esempio 2
Usare il metodo di Heun con h=0.2 per trovare i valori
approssimati della soluzione del problema ai valori
iniziali dell’esempio 1 in [0,1]. Confrontare gli errori con
quelli ottenuti tramite il metodo di Eulero. [excel]
Osservazione: gli errori sono notevolmente inferiori ad 1/10
di quelli ottenuti nell’esempio1(a). Ovviamente sono
necessari piu’ calcoli per ogni passo, ma il numero
richiesto di valutazioni di f(x,y) è solo doppio di quello
richiesto nell’esempio1(a). Come nel caso
dell’integrazione numerica, se f è complicata, è proprio il
numero di valutazioni della funzione a costituire gran
parte del ‘costo’computazionale del calcolo della
soluzione numerica.
98
Esercizio
Ricavare i valori in x=0.4,
x=0.8 e x=1 della
soluzione del problema ai
valori iniziali mediante i
metodi appena descritti
ed usando passi di
lunghezza h=0.2, h=0.1 e
h=0.05 per ciascun
metodo. Che cosa
indicano i risultati
riguardo ai valori della
soluzione in questi punti?
Confrontare i risultati con
la soluzione esatta
y=1/(1-x).
 y'  y

y
(
0
)

1

2
99
22 dic 09
Modello preda-predatore
(Lotka-Volterra)
 x' (t )  ax(t )  cx(t ) y (t )
 y ' (t )  by (t )  dx(t ) y (t )

a, b, c, d  0
100
Le equazioni di Lotka-Volterra, note anche come
equazioni preda-predatore, sono un sistema di
equazioni differenziali non lineari del primo
ordine. Tali equazioni forniscono un modello
matematico in grado di descrivere la dinamica di
un ecosistema in cui interagiscono soltanto due
specie animali: una delle due come predatore,
l'altra come la sua preda. Questa
modellizzazione matematica è stata proposta
indipendentemente da Alfred J. Lotka nel 1925 e
Vito Volterra nel 1926.
101
Supponiamo che nello stesso ambiente convivano
due sole specie di animali. Quelli della prima
specie sono erbivori, mentre quelli della
seconda, i predatori, sono carnivori e si nutrono
degli animali appartenenti alla prima specie, le
prede. Il nostro interesse è essenzialmente
rivolto allo studio dei legami tra le due specie:
come varia nel tempo il numero degli animali
appartenenti alle due specie? Possono le due
specie convivere nello stesso habitati oppure
una di esse è destinata ad estinguersi? Che
influenza hanno le caratteristiche biologiche
delle singole specie nell’equilibrio finale
dell’ecosistema?
102
Se centriamo la nostra analisi su questi problemi
possiamo, almeno in prima approssimazione,
trascurare alcuni aspetti biologici che sono
certamente connessi allo studio della evoluzione
delle due specie. Supponiamo, ad esempio, che
il cibo della prima specie sia disponibile sul
territorio senza limitazioni e ignoriamo ogni
differenziazione di età e di sesso tra gli animali.
E’ evidente che queste semplicifazioni portano a
considerare un ecosistema astratto. Nella
costruzione del modello viene quindi isolato
l’aspetto che piu’ interessa, in questo caso
l’interazione preda-predatore.
103
La validità del modello matematico come
modello di previsione dipende pertanto da
quanto i fattori trascurati siano ininfluenti
rispetto allo studio delle interazioni predapredatore.
Supponiamo che sia le prede che i predatori
si riproducano con tassi di natalità costanti
e si estinguano con tassi di mortalità
costanti.
104
Quindi, indicati con x(t) ed y(t)
rispettivamente il numero delle prede e dei
predatori presenti sul territorio al tempo t,
in assenza di predazione l’evoluzione nel
tempo delle due specie è descritta dalle
due equazioni differenziali:
x’(t)= a x(t)
y’(t)=-b y(t)
con a>0, b>0.
105
Si tratta di due equazioni, ciascuna a variabili
separabili, non collegate fra di loro. La soluzione è
x’(t)= x(0)*exp(a t)
y’(t)= y(0)* exp(-b t)
Il che vuol dire che il numero delle prede tende a
crescere indefinitamente, coerentemente con
l’ipotesi che il loro cibo è disponibile senza
limitazioni, mentre il numero dei predatori decresce
verso zero perche’ non si è tenuto conto della
attività di predazione, che è l’unica forma di
sostentamento della seconda specie.
106
Per descrivere l’attività di predazione
supponiamo che il numero delle prede
diminuisca di un fattore proporzionale a
x(t)*y(t)
Questo numero si interpreta come numero
di ‘incontri’ tra le due specie.
Supponiamo che il numero dei predatori
cresca in maniera analoga.
107
Si ottiene il seguente sistema:
 x' (t )  ax(t )  cx(t ) y (t )
 y ' (t )  by (t )  dx(t ) y (t )

a, b, c, d  0
La determinazione esplicita della soluzione di
questo sistema non è facile, a causa del legame
dato dal termine non lineare x(t)y(t).
108
E’ possibile effettuare un’analisi qualitativa
che fornisca in maniera semplice qualche
informazione sull’evoluzione delle due
specie.
Il primo passo consiste nel determinare i
punti di equilibrio, cioè quelli in cui la
derivata prima si annulla.
Omettiamo, per semplicità di scrittura, la
dipendenza da t e risolviamo il sistema
x’=0, y’=0.
109
Partiamo dal sistema:
x’(t)= a x(t) – c x(t)*y(t)
y’(t)= -b y(t) + d x(t)*y(t)
Imporre x’=0 ed y’=0 vuol dire chiedere
a x – c x*y=0
-b y + d x*y=0
Metto in evidenza:
x( a – c y) = 0
y(-b + d x) = 0
Verifico che E1=(0,0) è soluzione e che E2=(b/d,a/c) è
soluzione.
Ho quindi due punti di equilibrio.
110
Questo significa che se il numero iniziale delle
prede e dei predatori è dato dalle coordinate di
E1 (oppure dalle coordinate di E2), le due
popolazioni (preda e predatore) conservano nel
tempo esattamente lo stesso numero di membri.
Naturalmente questa situazione è assai
improbabile e ci interessa piuttosto sapere se,
partendo da una arbitraria situazione iniziale
(x(0), y(0)), l’evoluzione converga o meno ad uno
degli equilibri.
Analizziamo ora il comportamento qualitativo del
sistema per una generica scelta dei dati iniziali
x(0), y(0) strettamente positivi.
111
La situazione è descritta in figura, che rappresenta il segno delle componenti nel
campo vettoriale f1(x,y)=ax-cxy ; f2(x,y)=-by+dxy nel piano delle fasi. Se si parte da un
punto sull’asse delle x le prede tenderanno a moltiplicarsi perché non ci sono
predatori sul territorio. Se invece la condizione iniziale è sull’asse delle ordinate, allora
i predatori tenderanno ad estinguersi in assenza di prede e quindi la soluzione del
sistema convergerà all’equilibrio E1
112
E’ interessante osservare che non appena x o y sono uguali a zero anche le
corrispondenti derivate x’ e y’ si annullano;: se il numero dei membri di una delle
popolazioni è nullo ad un certo tempo, rimarrà nullo per tutti i tempi successivi.
Partendo invece da un punto con componenti strettamente positive, il segno delle
derivate x’ ed y’ varierà a seconda che ci si trovi nelle regioni I, II, III, o IV.
113
Quali tecniche conosciamo per studiare un
sistema di equazioni non lineare?
1. Si linearizza e si studia il sistema
linearizzato.
2. Tecniche alternative.
114
1. Dal sistema al sistema linearizzato
 x' (t )  ( a  cy) x  f1 ( x, y)

 y' (t )  (b  dx) y  f 2 ( x, y)
La stabilità dei punti di equilibrio può essere determinata
linearizzando il sistema utilizzando le derivate parziali.
La matrice jacobiana del modello preda-predatore è
 f1

x

J ( x, y ) 
 f 2
 x

f1 

 cx 
y   a  cy

 
f 2   dy
 b  dx 

y 
115
La matrice jacobiana
calcolata nel primo
punto di equilibrio è:
a 0 

J (0,0)  
 0  b
Gli autovalori sono uguali, rispettivamente, ad a e –b.
Siccome a>0 e b>0, si ha che gli autovalori i segni dei due
autovalori sono sempre diversi. Dunque il punto di equilibrio
nell'origine è un punto di sella. La stabilità di questo punto è
importante: se fosse stabile, valori di popolazione diversi da
zero potrebbero essere attratti da esso, e perciò la dinamica
del sistema porterebbe all'estinzione di entrambe le specie
per molti valori iniziali delle popolazioni. Visto però che il
punto è di sella, l'equilibrio è instabile, e l'estinzione di
entrambe le specie è quindi difficile (in effetti, può accadere
solo se le prede vengono estinte completamente in modo
artificiale, provocando la morte dei predatori a causa della
mancanza di cibo. Se invece sono i predatori ad essere
estinti, la popolazione delle prede cresce senza limite in
questo semplice modello).
116
La matrice jacobiana calcolata nel primo punto di
equilibrio è:
a
b  

c
  0
 cx   a  c c
 a  cy
d 
  
J ( x, y )  
b   ad
 b  dx   d a
 dy
b d  

c
d  c

bc 
 
d 
0 

Gli autovalori si ottengono annullando il determinante:
bc 

   
d   2  bc  ad
det 
d c
 ad   


c


Dato che gli unici valori che annullano il determinante sono
complessi, il punto è di equilibrio stabile. La parte reale è
uguale a zero in entrambi i casi, quindi il punto di equilibrio è
un centro. Ciò significa che il livello dei predatori e delle prede117
è ciclico, ed oscilla attorno a questo punto fisso.
2. Tecniche alternative
118
 x' (t )  ( a  cy) x  f1 ( x, y)

 y' (t )  (b  dx) y  f 2 ( x, y)
119
120
y
Proiezione sul
piano (x,y) delle
curve di livello
H(x,y)
121
x
122
Per analizzare il comportamento del sistema fuori dai punti di
equilibrio si può dimostrare che esso ammette la costante del
moto H(x,y). Dunque tutte le traiettorie del sistema nello spazio x
y giacciono sulle curve di livello della funzione H. La funzione ha
un minimo corrispondente al punto di equilibrio stabile ed è
convessa sul primo quadrante, da questo si può dedurre che in
una situazione generica con due popolazioni iniziali x e y il
sistema ha un comportamento oscillante che torna
periodicamente nello stato iniziale, con oscillazioni anche molto
grandi.
Le soluzioni non possiedono una espressione semplice in
termini delle funzioni trigonometriche. Se si approssimano le
soluzioni che si trovano in un intorno del punto di equilibrio
stabile linearizzando il sistema si ottiene un moto armonico
semplice in cui la popolazione dei predatori precede quella delle
prede con uno sfasamento uguale a π / 2.
123
124
Domande per l'esame
- Elementi propedeutici di matematica generale
o enunciato del teorema di Bolzano-Weierstrass
o teorema di Rolle
o teorema di Lagrange
o formula di Taylor; dimostrazione della convergenza del resto di ordine uno
o definizione di max, sup: quale relazione li lega?
o equazione di una retta per due punti
- Metodi per la risoluzione di equazioni differenziali
o Illustra il metodo di Eulero
o Mostra il modo in cui viene ricavato il metodo a partire dalla formula di Taylor
o Calcola l'errore ad un passo
o Calcola l'errore globale
o Dimostra l'ordine di grandezza dell'errore rispetto al passo.
o Scrivi il generico passo ricavato con il metodo di Eulero per la risoluzione dei
problemi di Caucy:
A. x'(t)=t^2; x(0)=0;
B. x'(t)=x(t); x(0)=1;
o Illustra il modo in cui si ricavano metodi di ordine superiore a partire dalla formula
di Taylor
o Illustra il modo in cui si ricava il metodo di Heun a partire dalla formula di Taylor
125
o Descrivi il generico passo della iterazione del metodo di Runge-Kutta ed il suo
ordine di grandezza (solo enunciato).
- Metodi per la risoluzione di integrali
o Illustra il metodo dei rettangoli
o Calcola l'errore del metodo dei rettangoli
o Dimostra l'ordine di grandezza dell'errore
o Illustra il metodo dei trapezi
o Calcola l'errore del metodo dei trapezi
o Dimostra l'ordine di grandezza dell'errore
o Illustra il metodo di Simpson
- Legame tra metodi per la risoluzione di equazioni differenziali e metodi per la
risoluzione di integrali
o Illustra la problematica (in generale)
o Mostra come il metodo di Eulero si riduce al metodo dei rettangoli nel caso di
problemi di Cauchy in cui
l'equazione differenziale è del tipo x'(t)=f(t)
o Mostra come il metodo di Heun si riduce al metodo dei trapezi nel caso di
problemi di Cauchy in cui
l'equazione differenziale è del tipo x'(t)=f(t)
- Modello preda-predatore
o E' un sistema di equazioni lineari?
o Quali tecniche sono disponibili per la sua risoluzione?
o Cos'e' un integrale primo?
o Quali particolarità hanno le traiettorie?
o Come si riduce il sistema in assenza di attività predatoria?
126