1

LUISSCorso di

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

Lezione del xx/11/2009Calcolo approssimato di integrali - 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

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).

Xxtx

ttftx

00 )(

)()(' T

0

0

)()( xdftxt

t

5

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

a

aFbFdttf )()()(

6

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)

7

Obbiettivo: calcolare UN VALORE APPROSSIMATO

dell’integrale

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).

b

a

dttf )(

8

Punto fondamentale

Calcolare a priori l’errore commesso con la specifica procedura adottata. Serve per garantire l’affidabilità del metodo e bisognerà pertanto stimare

b

a

b

a

dttgdttfE )()(

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

NkNabkatk ,,1,0,/)(

Il numero h=(b-a)/N è il passo della discretizzazione

9

METODO DEI RETTANGOLI

Devo calcolare l’area sottesa dal grafico della funzione.

a b

10

METODO DEI RETTANGOLI

Idea: costruisco la partizione di [a,b]

NkNabkatk ,,1,0,/)(

a=t0 t1 t2 b=tN

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

11

a=t0 t1 t2 b=tN

Considero i punti intermedi in ciascun intervallo

12

a=t0 t1 t2 b=tN

Considero il valore della funzione in ciascun punto intermedio in ciascun intervallo

13

a=t0 t1 t2 b=tN

Considero la funzione costante a tratti, avente, su ciascun intervallo, il valore costante uguale all’altezza della funzione calcolata nel punto intermedio dell’intervallo

Nkttttt

ftg kkkk ,,2,1),,[,

2)( 1

1

14

Considero la funzione costante a tratti, avente, su ciascun intervallo, il valore costante uguale all’altezza della funzione calcolata nel punto intermedio dell’intervallo

a=t0 t1 t2 b=tN

Nkttttt

ftg kkkk ,,2,1),,[,

2)( 1

1

Ovviamente g coincide con f in tutti i punti intermedi di ciascun intervallo.

15

a=t0 t1 t3 b=tN

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

abxxfI

N

k

kk

1

0

1

2

16

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).

17

Stima dell’errore

b

a

b

a

dttgdttfE )()(

18

MEMO: da matematica generale Proprietà dell’integrale definito

a b

f(x)

x

y

0

L’area si può calcolare dividendola in due parti: (a, c) e (c, b)

( ) ( ) ( )b c b

a a c

f x dx f x dx f x dx

c

19

Applico la proprietà dell’integrale ripetutamente:

1

0

1

1

2

1

1

)(

)()()()(

N

k

t

t

b

t

t

t

t

a

b

a

k

k

N

dttf

dttfdttfdttfdttf

Posso ripetere il procedimento per g e quindi

1

0

1

)()()()(N

k

t

t

b

a

b

a

k

k

dttgtfdttgdttfE

20

MEMO

baba In generale

dcbadcba

Il valore assoluto di una somma è minore od uguale alla somma dei valori assoluti

21

1

0

11

0

1

0

11

1

2)()()(

)()()()(

N

k

t

t

kkN

k

t

t

N

k

t

t

b

a

b

a

k

k

k

k

k

k

dttt

ftfdttgtf

dttgtfdttgdttfE

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

1

0

11

2)('

N

k

t

t

kkk

k

k

dttt

tf

Ho applicato la proprietà del valore assoluto

Scrivo l’espressione di g

22

Osservazione

),[,22 1

1

kkkk ttt

httt

tk tk+1(tk+tk+1)/2

t

Un punto t dista dal centro dell’intervallo meno di metà della lunghezza dell’intervallo.

23)(

2)('sup)(

2)('sup

)(2

)('sup)(2

)('sup

12

)('sup2

)('sup

2)('

2)('

1

01

1

0

1

0

1

0

1

0

1

11

11

abh

fNhh

f

Nhh

ftth

f

dth

fdth

f

dth

fdttt

tf

kba

kba

kba

N

kkkk

ba

N

k

t

t

kba

N

k

t

t

kba

N

k

t

t

k

N

k

t

t

kkk

k

k

k

k

k

k

k

k

),[,22 1

1

kkkk ttt

httt

Utilizzo questa proprietà

per svolgere il passaggio

24

2)()('sup)(

habfhE k

ba

Conclusione: ho dimostrato che

25

METODO DEI TRAPEZI

Nel metodo dei rettangoli si approssima f con g costante a tratti

a=t0 t1 t2 b=tN

26

METODO DEI TRAPEZI

Nel metodo dei trapezi si approssima f con g costruita partendo dalle rette secanti in ciascun intervallo

a=t0 t1 t2 b=tN

Quindi la funzione g è lineare in ogni intervallo e coincide con f negli estremi di ciascun intervallo.

27

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 b=tN

Quindi la funzione g è lineare in ogni intervallo e coincide con f negli estremi di ciascun intervallo.

28

METODO DEI TRAPEZIL’approssimazione corrisponde al calcolo dell’area evidenziata in

rosso.

Obbiettivi:

•calcolo dell’area

•stima dell’errore

29

Osservazione: rispetto alla normale visualizzazione, i

trapezi sono ruotati di 90 gradi:

MEMO: area di un trapezio

(base minore + base maggiore)*altezza /2

Base maggiore

Base minore

altezza

30

Equazione di una retta che passa per i punti

(x1,y1) e (x2,y2)

2 1

2 1

y ym

x x

2 1 1

2 1 1

y y y ym

x x x x

x1 x x2

y2

y

y1

x

y

da cui

1 2 1

1 2 1

y y y y

x x x x

da cui

2 11 1

2 1

( )y y

y y x xx x

da cui

2 11 1

2 1

( )y y

y x x yx x

MEMO: da matematica generale

31

Calcolo della retta secante nel singolo intervallo

2 11 1

2 1

( )y y

y x x yx x

ttk t tk+1

f(tk+1)

f(t)

f(tk)

……

Applico la formula

Considerando le particolari coordinate dei punti:

1,...,1,0),,[

)()()()(

)(

1

1

1

Nkttt

tftttt

tftftg

kk

kkkk

kk

Osservo che la funzione g coincide con la funzione f negli estremi dei singoli intervalli ed è lineare in ogni intervallo.

32

Area approssimata = somma delle aree dei singoli trapezi

2/)()(

2/)()()()(

1

01

1

011

N

kkk

k

N

kkkk

b

a

N

abtftf

tttftfdttf

Base minore+base maggiore altezza

Perché gli intervalli hanno tutti la stessa ampiezza

33

Stima dell’erroreIpotesi in più: f derivabile due volte in (a,b) con derivata f’’ limitata.

Calcolo l’errore nel singolo intervallo, poi sommo per ottenere l’errore totale.

Inizio il calcolo partendo dall’espressione di g

1,...,1,0),,[

)()()()(

)(

1

1

1

Nkttt

tftttt

tftftg

kk

kkkk

kk

Mfba

)(''sup),(

34

Inoltre g ed f coincidono sui punti della partizione

)()()(

)()()()(1

1k

kk

kkkk tt

tt

tftftftgtgtg

))((')()( kkk ttftftf

)()()(

)('))()(())()((1

1k

kk

kkkkk tt

tt

tftfftgtgtftf

Per il teorema di Lagrange

Sottraggo

Da cui

)()()(

)('))()(())()(()()(1

1k

kk

kkkkk tt

tt

tftfftgtgtftftgtf

Ora cerco di semplificare questa quantità

35

kk

kkk tt

tftff

1

1 )()()('

)()( 1 kk tftf

kkkkk ttftftf 11 )(')()(

kkkkkk

kkk ttff

tt

tftff

1

1

1 )(')(')()(

)('

Considero quindi

Lavoro dapprima su

Applico la formula di Lagrange ed ottengo

Sostituisco ed ottengo:

36

kkkkkk

kkk ttff

tt

tftff

1

1

1 )(')(')()(

)('

kkkk

kkkk

ttfttf

fff

11 )('')(''

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

Applico il teorema di Lagrange alla funzione derivata prima:

Sostituisco ed ottengo che in ciascun intervallo

2

),(

221

11

1

)(''sup)('')(''

)()()(')()(

hfhfttf

tttt

tftfftgtf

bakk

kkkk

kkk

37

Errore totale: sommo sugli N intervalli

31

0

31

01

2

1

0

21

0

21

0

)(

)()(

)()()()(

111

NMhMhttMh

dtMhdtMhdttgtf

dttgtfdttgdttfE

N

k

N

kkk

N

k

t

t

N

k

t

t

N

k

t

t

b

a

b

a

b

a

k

k

k

k

k

k

3)( NMhhEE Quindi la stima cercata è:

38

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 h0.

Il passo successivo riguarda l’approssimazione mediante una funzione quadratica.

Passi ulteriori possibili riguardano l’approssimazione della funzione assegnata mediante polinomi di ordine maggiore.

39

Interpolata quadratricaConsidero una funzione in [-1,1] e la sua interpolata

quadratica definita da:

]1,1[,)( 2 tCBtAttgDove le tre costanti sono determinate dalle tre condizioni

)1()1(

)0()0(

)1()1(

fCBAg

fCg

fCBAg

che chiedono che le funzioni coincidano su quei tre punti.

Osservo che l’unica soluzione del sistema lineare è:

)0(,2

)1()1(),0(

2

)1()1(fC

ffBf

ffA

40

Interpolata quadratricaQuindi

3

)1()0(4)1(

23

2]

23[

23

23)(

1

1

231

1

21

1

fff

CACBA

CBA

Ctt

Bt

AdtCBtAtdttg

Questo calcolo si generalizza facilmente ad un arbitrario intervallo [c,d] tramite un cambio di variabile

tcddc

t22

Che trasforma l’intervallo [-1,1] in [c,d]. L’integrale diventa quindi

d

c

dfdcfcfcdCAdttgdfdcfcf 6/)()2/)((49(23

2)()(2/)(4)(

41

Su un intervallo qualsiasi [c,d]

con semiampiezza dell’intervallo h h=(m-c)=(d-m)

)()(4)(3

)( dfmfcfh

dttgd

c

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”)

g(x)f(x)

h hc dm

42

Su più intervalli consecutiviApplico 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 b=tN

Quindi la funzione g è lineare in ogni intervallo e coincide con f negli estremi di ciascun intervallo.

43

Su più intervalli consecutiviApplico 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 b=tN

h h h h h h

2

021 )()(4)(

3

1)(

N

parikk

kkk

b

a

tftftfN

abdttf

(h=(b-a)/N)

44

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

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.

0

)( dttf