LUISS
Corso di
Metodi matematici per economia e finanza.
Prof. F. Gozzi
Lezione del xx/11/2009
Calcolo approssimato di integrali - 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.
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
4
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
5
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)
6
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).
7
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
8
METODO DEI RETTANGOLI
a
b
Devo calcolare l’area sottesa dal grafico della funzione.
9
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
10
Considero i punti intermedi in ciascun intervallo
a=t0
t1
t2
b=tN
11
Considero il valore della funzione in ciascun punto intermedio in ciascun
intervallo
a=t0
t1
t2
b=tN
12
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
13
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
14
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
 xk  xk 1  b  a
f

2

 N
a=t0
t1
t3
b=tN
15
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(=(b-a)/N).
16
Stima dell’errore
E
b
b
a
a
 f (t )dt   g (t )dt
17
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
18
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
19
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
20
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
21
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.
22
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
23
Conclusione: ho dimostrato che
h
E (h)  sup f ' ( k ) (b  a)
2
a  b
24
METODO DEI TRAPEZI
Nel metodo dei rettangoli si approssima f con g costante a tratti
a=t0
t1
t2
b=tN
25
METODO DEI TRAPEZI
Nel metodo dei trapezi si approssima f con g costruita partendo
dalle rette secanti in ciascun intervallo
a=t0
t1
t2
Quindi la funzione g è lineare in ogni intervallo e coincide con f negli
estremi di ciascun intervallo.
b=tN
26
METODO DEI TRAPEZI
Il metodo si chiama ‘dei trapezi’ perché approssimo l’area tramite
trapezi (eventualmente degeneri) invece che tramite rettangoli.
a=t0
t1
t2
Quindi la funzione g è lineare in ogni intervallo e coincide con f negli
estremi di ciascun intervallo.
b=tN
27
METODO DEI TRAPEZI
L’approssimazione corrisponde al calcolo dell’area evidenziata in
rosso.
Obbiettivi:
•calcolo dell’area
•stima dell’errore
28
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
29
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
30
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.
31
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
32
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
33
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à
34

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


35

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
36
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
2
k 0
 dt 
tk
  Mh 2 (t k 1  t k )  Mh 3 NMh 3
Quindi la stima cercata è:
E  E (h)  NMh
3
37
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.
38
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
39
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
f (c  d ) / 2  f (d ) g (t )dt  A  2C  d  c  f (c9  4 f ((c  d ) / 2)  f (40d ) / 6
3
c
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”)
41
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
42
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
43
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
0
Si può anche dimostrare che il metodo di Eulero modificato si
riduce a quello dei rettangoli, il metodo di Runge-Kutta a
quello di Simpson.
44