1

LUISSCorso di

Metodi matematici per economia e finanza. Prof. F. Gozzia.a. 2009/10

Integrazione numerica di equazioni differenziali- Prof.ssa G. Rotundo

2

Elenco degli argomentiFormule 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 integraliMetodo dei rettangoliMetodo dei trapeziMetodo di Simpson 1/3

Corrispondenza tra metodi di calcolo numerico di EDO ed integraliCorrispondenza tra il metodo di Eulero (modificato) ed il metodo dei rettangoliCorrispondenza tra il metodo di Heun ed il metodo dei trapezi

Modello preda-predatoreDescrizione al modelloFormalizzazione tramite un sistema di equazioni differenzialiLinearizzazioneIntegrale primo, grafici di curve di livello e traiettorie.

3

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

4

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.

16 nov 09

5

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

6

Il metodo di Eulero

Leonhard Euler (1707-1783)

7

Il metodo di EuleroE’ il metodo più semplice. Per illustrarlo si parte dalla definizione

di derivata:

)(')()(

lim0

txh

txhtxh

Per h abbastanza piccolo il rapporto incrementale è quasi uguale alla derivata:

))(,()()(

txtfh

txhtx

))(,()()( txtfhtxhtx da cui

)(')()(

txh

txhtx

))(,()(' txtftx Sostituendo

si ha che

8

))(,()()( txtfhtxhtx

))(,()()( txtfhtxhtx

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

La relazione

9

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

10

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

t0 t1t2

x0

x1

x2

La funzione approssimante, lineare a tratti, può essere descritta, t: 0tT, dalla espressione: ...2,1,0),)1(,[),,()()( khkkhtxtfhktxtx kkh

11

Esempio 1

Si verifica che soddisfa una equazione differenziale, infatti:

2)( ttx Considero la seguente funzione:

ttxtx

txtxt

ttxt

ttx

/)(2)('

)(2)('

2)('

2)('2

ttxtxtf

txtftx

/)(2))(,(

))(,()('

2)( ttx Quindi è soluzione della equazione differenziale:

con

12

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 è: ))(,()()( txtfhtxhtx

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: x’(t)=f(t,x(t))))(,()()( txtfhtxhtx

kkkk txxx /22.01

13

Confronto con lo 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

da cui ricavo una stima dell’errore dell’approssimazione di x(h) tramite x(0)+hx’(0) :

Ammesso che esista

(*)

14

Esempio 2

Considero la seguente funzione:

Si verifica che soddisfa

un problema di Cauchy: )()('

)('

)(

txtx

etx

etxt

t

1)0(

)()('

x

txtx

15

Costruzione della approssimazione numerica della soluzione di

1)0(

))(()()('

x

txftxtx

hxhxhxhxx 1)0(1)0(')0()(1

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

Scrivo i punti successivi che si ottengono applicando il metodo:

Osservo che la variabile t non appare da sola, ma soltanto tramite x, quindi l’approssimazione di Eulero diventa:

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

Scrivo il primo punto che si ottiene tramite l’approssimazione con il metodo di Eulero:

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

Quindi l’approssimazione di Eulero diventa: kkk xhxx 1

16

Quindi la funzione soluzione del problema è:tetx )(

khx k )1( e le iterate approssimanti sono:

0 0.2 0.4 0.6 0.8 1 1.2 1.4 1.61

1.5

2

2.5

3

3.5

4

4.5

5

t

x(t)

tetx )(

0x1x

2x

3x

4x

k t x(t) xk

0 0.00 1.00 1.00

1 0.40 1.49 1.40

2 0.80 2.23 1.96

3 1.20 3.32 2.74

4 1.60 4.95 3.84

Osservo che anche in questo caso la distanza (barre verticali in nero) della soluzione da quella approssimata aumenta al crescere di k.

17

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

Programma matlab eulero.m

18

Nota

Xxtx

ttxftx

00 )(

))(()(' T

Da questo momento in poi consideriamo, per semplicità, il problema autonomo

19

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

hxhxx

xx

kkk

)(),( kk xfhx Il metodo di Eulero corrisponde ad una particolare scelta di :

20

Analisi della qualità del metodo

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

2. l’errore globale di approssimazione (devo esaminare cosa succede se itero

l’approssimazione N volte ed N∞)

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

21

1. (a)Misura della qualità di un metodo ad un passo

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

),(),(),( 000 hxhxhxe

0se)),0((

0se,)0()(

),( 0

hxf

hh

xhxhx

con

22

1.(b) Dimostro che 0)),0((lim0

hxeh

Svolgimento:

0lim

))(()0()(lim

))(()0()(lim

))(()0()(

lim)),0((lim

2

00

0

00

h

Ch

h

txhfxhx

h

txhfxhx

txfh

xhxhxe

hh

h

hh

utilizzo (*)

23

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.

24

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.

25

La funzione soluzione del problema è:tetx )(

khx k )1( e le iterate della approssimazione di Eulero sono:

Pongo t=hk, così stimo l’errore solo nei punti della successione tk .

)1log(

)1()()(hkkh

kkhkkk

ee

hexkhxxtx

Voglio stimare l’errore, cioè la differenza tra la soluzione esatta e quella approssimata.

Devo semplificare questa espressione in modo da capirne facilmente l’andamento quando il numero k di passi aumenta (all’infinito).

1)0(

)()('

x

txtx

2.(a) Illustrazione del problema tramite un esempio

Considero nuovamente il problema di Cauchy

26

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

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.

27

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:

28

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.

29

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[

Definizione: l’errore globale della approssimazione è dato da:

30

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

23 nov 09

31

;)('sup;)(sup tfLtfMRtRt

2/)1()( LTN eMhTxx

dhkxfhkxkhxh

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

Dimostrazione

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

stt

s

tsdsxfsxdxfsxtx0

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

)(khxxd kk Ponendo si ha dunque

dhkxfhkx

xhfxdh

kkk

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

)(

0

11

Caso particolare:

32

h

kk

h

kkk

dhkxfxfd

dhkxfhkx

xhfxd

0

11

0

11

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

ora sommo e sottraggo la stessa quantità

33

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

stt

s

tsdsxfsxdxfsxtx0

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

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

34

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

q.e.d

35

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

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)

36

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

37

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

38

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

)(!3

)0("'

2

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

32 hRx

hx

hxhxhx

),(),(),( 000 hxhxhxe

0se)),0((

0se,)0()(

),( 0

hxf

hh

xhxhx

Ricordo anche l’espansione in serie di Taylor

Da cui

)(!3

)0("'

2

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

32 hRx

hx

hxhxhx

Sostituendo l’espressione di (x0,h) :

),()(

!3

)0("'

2

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

h

hRxh

xhxhte

h

hRhxe

)(),( 3

0 )(),( 30 hOhxe

quindi

39

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.

40

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.

41

Ipotesi: derivabile due volte

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

)(2

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

2 hRx

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

42

Esempio

))(()(),( 000 xChfxBfxAfhx

0)0,(')0,(2

)()(')( 00

000 xhx

xfxfhxf n

)()(')()(2

)()(')( 00

000 xfxhBCfxfBA

xfxfhxf

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

43

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/2Corrisponde ad utilizzare come pendenza la pendenza calcolata nel punto

intermedio di ciascun intervallo della partizione.

44

Metodo di Runge-Kutta

Si ottiene definendo

6/)22(),( 43210 kkkkhx

Con

)(

)2/(

)2/(

)(

304

203

102

01

hkxfk

hkxfk

hkxfk

xfk

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

45

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 modificatoE(h)=O(h4) per il metodo di Runge-Kutta

46

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

47

Metodi numerici per il calcolo approssimato di integrali

30 nov 09

•Metodo dei rettangoli

•Metodo dei trapezi

•Metodo di Simpson

48

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

49

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

50

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

51

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

52

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

53

METODO DEI RETTANGOLI

Devo calcolare l’area sottesa dal grafico della funzione.

a b

54

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

55

a=t0 t1 t2 b=tN

Considero i punti intermedi in ciascun intervallo

56

a=t0 t1 t2 b=tN

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

57

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

58

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.

59

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

abttfI

N

k

kk

1

0

1

2

60

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

61

Stima dell’errore

b

a

b

a

dttgdttfE )()(

62

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

63

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

64

MEMO

baba In generale

dcbadcba

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

65

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

66

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.

67)(

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

68

2)()('sup)(

habfhE k

ba

Conclusione: ho dimostrato che

69

METODO DEI TRAPEZI (o di Crank-Nicolson)Mentre nel metodo dei rettangoli

si approssima f con g costante a tratti,

a=t0

t1 t2 b=tN

a=t0

t1 t2 b=t

N

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

Il metodo si chiama ‘dei trapezi’ perché approssimo l’area tramite trapezi (eventualmente degeneri) invece che tramite rettangoli.

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

70

L’approssimazione corrisponde al calcolo dell’area evidenziata in rosso.

Obbiettivi:

•calcolo dell’area

•stima dell’errore

Osservazione: alcuni trapezi possono essere degeneri, per esempio

composti da triangoli o rettangoli

71

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

72

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

73

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.

74

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

75

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),(

76

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à

77

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:

78

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

79

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

)()( 23 abMhNMhhEE Quindi la stima cercata è:

80

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.

81

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

82

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

dfdcfcfcdCAdttg 6/)()2/)((4)(23

2)(

83

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

84

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.

85

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)

86

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

87

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

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.

0

)( dttf

88

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.

14 dic 09

89

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 integraliMetodo dei rettangoli

Metodo dei trapezi

Metodo di Simpson 1/3

Xxtx

ttxtftx

00 )(

))(,()(' T

Xxtx

ttftx

00 )(

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

Confronto tra metodi per la risoluzione di EDO e metodi per il calcolo di integrali

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.

901

2

1

0

kt

t

t

t

),()()(

),()()(

),()()(

)(

11

111122

000011

00

kkkkkk xtfhxhtxtxx

xtfhxhtxtxx

xtfhxhtxtxx

txx

Ricordo che il metodo di Eulero individua la successione approssimante è definita per ricorrenza da

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.

91

),()()( 11 kkkkkk xtfhxhtxtxx

2

002211000

2211000

221112223

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

),(),(),(

),()),((),(

lll xtfhxxtfxtfxtfhx

xtfhxtfhxtfhx

xtfhxtfhxxtfhxx

1

00

1100011000

110001112

0001

00

),(

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

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

)(

lll xtfhx

xtfxtfhxxtfhxtfhx

xtfhxtfhxxtfhxxxtfhxx

txx

k

lllkkkk xtfhxxtfhxx

001 ),(),(

scrivo xk+1 raggruppando i termini in modo diverso

(**) quindi

92

),()()( 11 kkkkkk xtfhxhtxtxx

1

0

),(N

lllN xtfhx

1

00 ),(

N

lllN xtfhxx

Quindi la trasformazione della approssimazione che si ottiene con il metodo di Eulero

Considero, per semplicità, il caso in cui x0=0 :

permette di scrivere

1

0

1

2

N

k

kk ttfhI

Ricordo l’approssimazione dell’integrale che si ottiene dal metodo dei rettangoli:

21kk

l

tttLe due formule coincidono se , e questo di fatto

significa che stiamo usando il metodo di Eulero modificato.

93

),(1 hxhxx kkk ))(()(),( kkkk xChfxBfxAfhx

)),(2

,2

(),,( kkkkkk xtfh

xh

tfhxt

Infatti, la modifica che viene considerata nel metodo di Eulero modificato corrisponde a scrivere il generico passo della iterazione

ponendo A=0, B=1, C=1/2, quindi ))(2

(),( kkk xfh

xfhx

)2

(),,(h

tfhxt kkk

k

lllkkkk hxthxhxthxx

001 ),,(),,(

Avendo posto, per semplicità, x0=0 , fissando k=N e sostituendo si ottiene:

con

Reintroducendo esplicitamente la dipendenza dal tempo si ha:

Poiché ora stiamo trattando il caso di funzioni che dipendono soltanto dal tempo l’espressione si semplifica e si ha:

Questo significa che xk+1 puo’ essere scritto, analogamente a quanto svolto in (**), come

1

0

)2

(N

llN

htfhx Ovvero la stessa approssimazione del metodo dei rettangoli,

che considera come altezza del rettangolo il valore dellafunzione nel punto intermedio di ciascun intervallo della partizione. Senza questa modifica, la formula di Eulero considera rettangoli con altezza data dal valore della funzione nel primo estremo di ciascun intervallo.

94

Corrispondenza tra il metodo di Heun ed il metodo dei trapezi

))(1(2

1)(

2

1),( kkkk xhfxfxfhx

),(1 hxhxx kkk

))(()(),( kkkk xChfxBfxAfhx

e quindi diventa

))(,(2

1),(

2

1),,( kkkkkkk xhfxhtfxtfhxt

2

)()()(

2

1)(

2

1),,(

htftfhtftfhxt kk

kkkk

ovvero il valor medio di f rispetto agli estremi del singolo intervallo della partizione.

con con A=1/2, B=1/2, C=1

Reintroducendo esplicitamente la dipendenza dal tempo si ha:

Poiché ora stiamo trattando il caso di funzioni che dipendono soltanto dal tempo l’espressione si semplifica e si ha:

Nel metodo di Heun l’iterazione è definita da:

95

Avendo posto, per semplicità, x0=0 , fissando k=N e sostituendo si ottiene:

1

0 2

)()(N

l

llN

htftfhx

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.

Ricordo il valore che si ottiene con il metodo dei trapezi:

1

0

11

01

1

011

2

)()(2/)()(

2/)()()()(

N

k

kkN

kkk

k

N

kkkk

b

a

htftf

N

abtftf

tttftfdttf

t t+h

Base minore+base maggiore altezza

Cioè esattamente ciò che si ottiene con il metodo di Heun

96

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

1)0(y

yxdx

dy

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)

21 dic 09

97

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

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

98

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.

99

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

1)0(

' 2

y

yy

100

Modello preda-predatore (Lotka-Volterra)

0,,,

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

dcba

tytdxtbytytytcxtaxtx

22 dic 09

101

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.

102

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?

103

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.

104

La validità del modello matematico come modello di previsione dipende pertanto da quanto i fattori trascurati siano ininfluenti rispetto allo studio delle interazioni preda-predatore.

Supponiamo che sia le prede che i predatori si riproducano con tassi di natalità costanti e si estinguano con tassi di mortalità costanti.

105

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.

106

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.

107

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.

108

Si ottiene il seguente sistema:

La determinazione esplicita della soluzione di questo sistema non è facile, a causa del legame dato dal termine non lineare x(t)y(t).

0,,,

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

dcba

tytdxtbytytytcxtaxtx

109

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.

110

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=0Metto in evidenza: x( a – c y) = 0

y(-b + d x) = 0Verifico che E1=(0,0) è soluzione e che E2=(b/d,a/c) è

soluzione. Ho quindi due punti di equilibrio.

111

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.

112

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

113

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.

114

Quali tecniche conosciamo per studiare un sistema di equazioni non lineare?

1. Si linearizza e si studia il sistema linearizzato.

2. Tecniche alternative.

115

La stabilità dei punti di equilibrio può essere determinata linearizzando il sistema utilizzando le derivate parziali.

La matrice jacobiana del modello preda-predatore è

dxbdy

cxcya

y

f

x

fy

f

x

f

yxJ22

11

),(

),()()('

),()()('

2

1

yxfydxbty

yxfxcyatx

1. Dal sistema al sistema linearizzato

116

La matrice jacobiana calcolata nel primo punto di equilibrio è:

b

aJ

0

0)0,0(

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

117

La matrice jacobiana calcolata nel primo punto di equilibrio è:

0

0),(

c

add

bc

d

bdb

c

ad

d

bc

c

aca

dxbdy

cxcyayxJ

Gli autovalori si ottengono annullando il determinante:

c

ad

d

bc

c

add

bc

2det

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 prede è ciclico, ed oscilla attorno a questo punto fisso.

118

2. Tecniche alternative

119

),()()('

),()()('

2

1

yxfydxbty

yxfxcyatx

120

121

x

yProiezione sul piano (x,y) delle curve di livello

H(x,y)

122

123

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.

124

125

Domande per l'esame

- Elementi propedeutici di matematica generaleo enunciato del teorema di Bolzano-Weierstrasso teorema di Rolle o teorema di Lagrange o formula di Taylor; dimostrazione della convergenza del resto di ordine unoo definizione di max, sup: quale relazione li lega?o equazione di una retta per due punti

- Metodi per la risoluzione di equazioni differenzialio Illustra il metodo di Euleroo Mostra il modo in cui viene ricavato il metodo a partire dalla formula di Tayloro Calcola l'errore ad un passoo Calcola l'errore globaleo 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 Tayloro Descrivi il generico passo della iterazione del metodo di Runge-Kutta ed il suo ordine di

grandezza (solo enunciato).

126

- Metodi per la risoluzione di integralio Illustra il metodo dei rettangolio Calcola l'errore del metodo dei rettangolio Dimostra l'ordine di grandezza dell'erroreo Illustra il metodo dei trapezio Calcola l'errore del metodo dei trapezio Dimostra l'ordine di grandezza dell'erroreo Illustra il metodo di Simpson

- Legame tra metodi per la risoluzione di equazioni differenziali e metodi per la risoluzione di integralio 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-predatoreo 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?