Integrazione equazioni Lotka Volterra con Runge-Kutta di ordine quattro

Michele Mancusi

Introduzione all’equazioni di Lotka-Volterra:

N.B. In fondo al documento potete trovare il istato del programma

È stato sviluppato il programma in questione al fine di studiare le Equazioni di Lotka-Volterra, note anche

come equazioni preda-predatore, le quali 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 Lotka nel 1925 e Volterra nel 1926.

Le equazioni in questione sono:

{

𝑑𝑥

𝑑𝑡= (𝑎 − 𝑏𝑦)𝑥

𝑑𝑦

𝑑𝑡= (−𝑐 + 𝑑𝑥)𝑦

𝑥(𝑡) rappresenta il numero di prede al tempo t, mentre 𝑦(𝑡) rappresenta il numero di predatori al tempo t.

Le derivate 𝑑𝑥

𝑑𝑡 e

𝑑𝑦

𝑑𝑡 rappresentano i tassi di crescita nel tempo delle popolazioni di prede e predatori, quindi

le velocità con cui variano del tempo le due popolazioni. I parametri a,b,c,d sono coefficienti positivi che

descrivono come interagiscono tra loro le due specie: a rappresenta il tasso di crescita naturale della

popolazione preda; b il coefficiente di impatto della predazione sulle prede. Per quanto riguarda i predatori,

invece, assumiamo che c rappresenti il tasso di mortalità dei predatori e d il coefficiente di efficienza della

predazione sulle prede.

Il modello di Lotka-Volterra formula una serie di ipotesi riguardanti l'ambiente e l'evoluzione delle

popolazioni di predatori e prede:

La popolazione preda trova ampio cibo a tutte le ore.

L' approvvigionamento alimentare della popolazione predatrice dipende interamente dalla dimensione della

popolazione predata.

La velocità di variazione della popolazione è proporzionale alle sue dimensioni.

Durante il processo, l'ambiente non cambia in favore di una specie e l'adattamento genetico è

sufficientemente basso.

I predatori hanno appetito senza limiti.

Così come sono utilizzate le equazioni differenziali, la soluzione è deterministica e continua. Questo, a sua

volta, implica che le generazioni sia dei predatori che delle prede sono continuamente sovrapposte.

Oltre al modello sopra descritto sono state studiate anche altre sue varianti con altrettanti sistemi di equazioni

differenziali simili al precedente.

Descrizione del programma:

Il programma integra il sistema di equazioni differenziali sopra descritto. Una volta avviato chiede all’

utente quale dei quattro modelli proposti desidera studiare: Modello classico, con piccole perturbazioni, con

competizione intraspecie, o modello con pesca indiscriminata. In seguito il programma chiede in che modo

l’utente desidera procede, cioè se i coefficienti sopra descritti, le condizioni iniziali e il passo di integrazione

debbano essere scelti dall’utente, utilizzati quelli predefiniti oppure scelti a caso in un intervallo compreso

tra 0 e 20. Dopo aver scelto a configurazione desiderata il programma stampa su schermo i punti equilibrio

del sistema scelto con le relative matrici jacobiane calcolate in tali punti.

Inoltre crea diversi file: uno (LV%.dat dove % può essere ‘std’,’sp’o ‘fish’) contiene i valori delle prede e

dei predatori ad ogni istante di tempo e l’hamiltoniana del sistema (costante del moto), un file (LVerr.dat) i

dati necessari allo studio dell’errore dell’algoritmo d’integrazione e un altro (LVfun.dat) i dati necessari per

poter graficare la funzione H(x,y) hamiltoniana del sistema. (Si avverte che a seconda della configurazione

scelta alcuni file o alcuni dati possono non essere prodotti).

I valori predefiniti sonno stati scelti in modo da far visualizzare, in maniera indicativa al lettore, una delle

possibilità in cui il sistema può evolvere, ma naturalmente non le esaurisce tutte quante.

Per studiare tale sistema è stato utilizzato come algoritmo di integrazione Runge-Kutta di ordine quattro. Ci

si aspetta dunque che l’errore commesso sulla quantità conservata H ad un determinato tempo t definito

come 𝐻(𝑡)−𝐻(𝑡0)

𝐻(𝑡0) in funzione di dt sia di grado 4, e dunque una retta con coefficiente angolare 4 su carta

doppio logaritmica.

L’algoritmo utilizzato è il seguente:

�̇� = 𝑓(𝑥, 𝑦, 𝑡)

�̇� = 𝑔(𝑥, 𝑦, 𝑡)

Dove

�̇� = 𝑥(𝑎 − 𝑏𝑦)

�̇� = 𝑦(−𝑐 + 𝑑𝑥)

𝑘1 = 𝑓(𝑥𝑛, 𝑦𝑛 , 𝑡𝑛),

𝑙1 = 𝑔(𝑥𝑛, 𝑦𝑛 , 𝑡𝑛),

𝑘2 = 𝑓(𝑥𝑛 +1

2𝑑𝑡 ∙ 𝑘1, 𝑦𝑛 +

1

2𝑑𝑡 ∙ 𝑙1, 𝑡𝑛 +

1

2𝑑𝑡),

𝑙2 = 𝑔(𝑥𝑛 +1

2𝑑𝑡 ∙ 𝑘1, 𝑦𝑛 +

1

2 𝑑𝑡 ∙ 𝑙1, 𝑡𝑛 +

1

2𝑑𝑡),

𝑘3 = 𝑓(𝑥𝑛 +1

2𝑑𝑡 ∙ 𝑘2, 𝑦𝑛 +

1

2𝑑𝑡 ∙ 𝑙2, 𝑡𝑛 +

1

2𝑑𝑡),

𝑙3 = 𝑔(𝑥𝑛 +1

2𝑑𝑡 ∙ 𝑘2, 𝑦𝑛 +

1

2 𝑑𝑡 ∙ 𝑙2, 𝑡𝑛 +

1

2𝑑𝑡),

𝑘4 = 𝑓(𝑥𝑛 + 𝑑𝑡 ∙ 𝑘3, 𝑦𝑛 + 𝑑𝑡 ∙ 𝑙3, 𝑡𝑛 + 𝑑𝑡),

𝑙4 = 𝑔(𝑥𝑛 + 𝑑𝑡 ∙ 𝑘3, 𝑦𝑛 + 𝑑𝑡 ∙ 𝑙3, 𝑡𝑛 + 𝑑𝑡),

𝑡𝑛+1 = 𝑡𝑛 + 𝑑𝑡

𝑥𝑛+1 = 𝑥𝑛 +1

6𝑑𝑡(𝑘1 + 2 ∙ 𝑘2 + 2 ∙ 𝑘3 + 𝑘4)

𝑦𝑛+1 = 𝑦𝑛 +1

6𝑑𝑡(𝑙1 + 2 ∙ 𝑙2 + 2 ∙ 𝑙3 + 𝑙4)

MODELLO CLASSICO.

Il primo modello studiato è il sistema di equazioni di Lotka-Volterra sopra descritto.

Si noti che in assenza dei predatori (y0=0), le prede crescerebbero esponenzialmente secondo il modello

maltusiano, infatti essendo y = 0 ⇒ �̇� = 𝑎𝑥 la cui soluzione è 𝑥 = 𝑥0𝑒𝑎𝑡, infatti integrando con il

programma sviluppato e ponendo come condizione iniziale y0=0 il risultato è una crescita esponenziale delle

prede.

Lo stesso vale per i predatori, infatti in assenza di prede essi si estinguerebbero in maniera esponenziale:

infatti ponendo x0=0 la soluzione dell’equazione riguardante i predatori risulta essere 𝑦 = 𝑦0𝑒−𝑐𝑡, indicando

un estinzione esponenziale, così come si può anche evincere dal grafico sotto riportato

.

Il sistema sopra descritto ha due punti critici. Il primo è nell’origine, mentre il secondo si trova nel primo

quadrante. Si ha un punto di equilibrio nelle popolazioni quando il livello di entrambe rimane costante, cioè

quando 𝑥(𝑡) = 𝑥0 e 𝑦(𝑡) = 𝑦0 per ogni tempo t, e di conseguenza le derivate 𝑥 ′(𝑡) e 𝑦 ′(𝑡) sono uguali a

zero. Sostituendo questi valori nelle equazioni si ottengono le seguenti relazioni: Per trovare tali punti

bisogna risolvere il sistema {0 = (𝑎 − 𝑏𝑦)𝑥

0 = (−𝑐 + 𝑑𝑥)𝑦 Risolvendo questo sistema di equazioni si trovano due

soluzioni per la coppia (𝑥0, 𝑦0) corrispondenti a due punti di equilibrio:

(𝑥0, 𝑦0) = (0,0) e (𝑥0, 𝑦0) = (𝑐

𝑑,

𝑎

𝑏)

Nel primo caso naturalmente se le due popolazioni hanno 0 individui allora continueranno ad avere 0

individui in ogni istante successivo. Il secondo corrisponde invece alla situazione in cui i predatori

incontrano e mangiano, in ogni unità di tempo, un numero di prede esattamente uguale al numero di prede

che nascono, e questo numero di prede corrisponde proprio alla soglia critica di cibo che fa rimanere

stazionaria la popolazione dei predatori.

La stabilità dei punti di equilibrio può essere determinata linearizzando il sistema utilizzando la matrice

jacobiana del modello preda-predatore:

𝐽(𝑥, 𝑦) = [𝑎 − 𝑏𝑦 −𝑏𝑥

𝑑𝑦 𝑑𝑥 − 𝑐]

La matrice Jacobiana calcolata nel primo punto di equilibrio (0,0) è:

𝐽(0,0) = [𝑎 00 −𝑐

]

Gli autovalori di questa matrice sono 𝜆1 = 𝑎 e 𝜆1 = −𝑐 . Dato che 𝑎 e 𝑐 sono quantità positive, i segni dei

due autovalori sono sempre diversi. Dunque il punto di equilibrio nell'origine è un punto di sella. Se fosse un

punto di equilibrio 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. Dal momento che il punto è di sella, l'equilibrio è instabile e l'estinzione di entrambe le specie è

quindi difficile (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).

La matrice jacobiana calcolata nel secondo punto di equilibrio è:

𝐽(𝑐

𝑑,𝑎

𝑏) = [

0 −𝑏𝑐

𝑑𝑎𝑑

𝑏0

]

Gli autovalori di questa matrice sono:

𝜆1 = 𝑖√𝑎𝑐 e 𝜆1 = −𝑖√𝑎𝑐

poiché la parte reale degli autovalori è nulla siamo in un caso particolare e non si possono trarre conclusioni

dall'analisi lineare. Tuttavia, il sistema ammette una costante del moto H(x,y) e le curve di livello di questa

funzione sono traiettorie chiuse che circondano il punto fisso. Di conseguenza, le traiettorie delle popolazioni

di prede e predatori, oscillano intorno a questo punto fisso. La funzione H ha un minimo corrispondente al

punto di equilibrio stabile ed è convessa sul primo quadrante, da cui si può dedurre che in una situazione

generica con due popolazioni iniziali 𝑥 e 𝑦 il sistema ha un comportamento oscillante che

torna periodicamente nello stato iniziale, con oscillazioni anche molto grandi. Da quanto appena detto

possiamo quindi concludere che il punto di equilibrio (𝑥0, 𝑦0) = (𝑐

𝑑,

𝑎

𝑏) risulta essere stabile.

Le soluzioni di questo sistema non possiedono un’espressione semplice in termini di funzioni

trigonometriche, ma si può osservare che nel tempo sia la popolazione delle prede che quella dei predatori

seguono un evoluzione periodica

Qui sopra è riportato il grafico dell’evoluzioni nel tempo delle popolazioni di prede e di predatori.

Nelle situazioni di non-equilibrio, dunque, si ha che i predatori prosperano quando c'è abbondanza di prede

ma, alla lunga, si ritrovano senza cibo sufficiente per sfamare l'intera popolazione e cominciano ad

estinguersi. Mentre la popolazione dei predatori decresce quella delle prede aumenta di nuovo; questa

dinamica continua in un ciclo di crescita e decrescita.

Si possono anche plottare le soluzioni senza rappresentare il tempo, ma con un asse che rappresenta il

numero delle prede e l'altro che rappresenta il numero di predatori. Le soluzioni sono curve chiuse, e vi è una

quantità che si conserva ad ogni curva:

𝐻(𝑥, 𝑦) = 𝑑𝑥 − 𝑐 ∙ 𝑙𝑛(𝑥) + 𝑏𝑦 − 𝑎 ∙ ln (𝑦)

Il grafico (x,y) sopra riportato mostra alcune orbite del modello di Lotka - Volterra corrispondenti a curve di

livello della funzione H(x,y), tutte le traiettorie del sistema nello spazio x-y giacciono sulle curve di

livello della funzione H.

Inoltre tramite il programma è stato possibile ricostruire la funzione H(x,y) in maniera tridimensionale

rappresentando sugli assi x ed y il numero di prede e predatori e sull’asse z il valore della costante H

corrispondente. Per fare ciò si sono lasciati invariati i coefficienti a,b,c,d,e si sono variate le condizioni

iniziali x0,y0, ricostruendo così le varie curve di livello e calcolando per ogni curva il valore della costante H

corrispondente.

Di seguito è riportato il grafico tridimensionale della funzione H(x,y).

COMPETIZIONE INTRASPECIE

In seguito è stato studiato un sistema leggermente diverso dal sistema Lotka-Volterra sopra descritto. Infatti

sono state introdotte delle piccole perturbazioni all’interno del sistema ed è stato visto come esso evolveva.

Il sistema così generato può essere espresso dalle seguenti equazioni differenziali:

{

𝑑𝑥

𝑑𝑡= 𝑥(𝑎 − 𝑏𝑦 − 𝜀1𝑥)

𝑑𝑦

𝑑𝑡= 𝑦(−𝑐 + 𝑑𝑥 − 𝜀2𝑦)

Questo sistema si presenta uguale al precedente ma con l’aggiunta, in entrambe le equazioni, di un termine:

rispettivamente −𝜀1𝑥 per la prima equazione −𝜀2𝑦 per la seconda. I parametri 𝜀2 e 𝜀1 possono essere

considerati come dei coefficienti di interferenza nelle due specie, i quali interferiscono con il normale

andamento della popolazione determinando una progressiva riduzione della popolazione.

In questo caso i parametri 𝜀2 e 𝜀1 rappresentano dei coefficienti di competizione all’interno della stessa

specie: infatti mentre ci sono sempre due popolazioni, una di prede l’altra di predatori, che interagiscono tra

loro normalmente (una cibandosi dell’altra), esiste all’interno della stesa specie una concorrenza che

danneggia gli individui della stessa popolazione, quindi la specie stessa e questi due coefficienti

rappresentato quanto feroce possa essere la concorrenza, sia all’interno delle prede, sia all’interno dei

predatori, entrambi ad esempio nel procurarsi il cibo, o a sopravvivere. Inoltre nel sistema classico la

popolazione, come abbiamo già detto, crescerebbe esponenzialmente in assenza dei predatori. In questo

nuovo sistema, la crescita esponenziale è stata sostituita con una più verosimile crescita logistica(se il

termine 𝜀1 non risulta essere troppo grande): All'inizio la crescita è quasi esponenziale, successivamente

rallenta, diventando quasi lineare, per raggiungere una posizione asintotica dove non c'è più crescita.

Anche in questo modello i coefficienti a,b,c,d hanno lo stesso significato biologico che avevano nel modello

classico Lotka-Volterra.

In questo caso non è nata una quantità conservata o costante del moto.

Di seguito è riportato il grafico del numero di prede in funzione del tempo, esso rappresenta una crescita

logistica delle prede in assenza di predatori (y0=0).

Di seguito è riportato il grafico dei predatori in funzione delle prede sviluppato a partire dal sistema

lievemente perturbato sopra riportato (ponendo 𝜀2 e 𝜀1 pari a 0.0001). Come si può notare la figura risulta

avere un contorno più “spesso”: infatti queste pertubazioni, seppur lievi, tendono a far tendere il sistema

verso un punto di equilibio, anche se molto lentamente, traciando quindi una leggera spirale.

I punti di equilibrio sono stati ottenuti come nel modello precedente ponendo entrambe le derivate uguali a

zero, i punti di equilibrio ottenuti sono:

(𝑥, 𝑦) = (𝑎𝜀2+𝑏𝑐

𝑏𝑑+𝜀2𝜀1,

𝑎𝑑−𝑐𝜀1

𝑏𝑑+𝜀2𝜀1) se solo se 𝑏𝑑 + 𝜀2𝜀1 ≠ 0 𝑒 𝑏 ≠ 0;

(𝑥, 𝑦) = (0,0);

(𝑥, 𝑦) = (0, −𝑐

𝜀2 ) ;

(𝑥, 𝑦) = (𝑎

𝜀1, 0).

Senza ricorrere allo studio della matrice jacobiana cerchiamo di vedere se questo è un punto di

equilibrio stabile:

Dal grafico possiamo vedere che la traiettoria è ciclica intorno al primo punto di equilibrio, verso cui

converge, questo significa che il punto di equilibrio è stabile (anche senza studiare la matrice jcobiana).

Questa seconda configurazione è stata ottenuta aumentato il valore delle perturbazioni (ponendo 𝜀2 e 𝜀1 pari

a 0.01) ottenendo così una figura a spirale molto più pronunciata della precedente, in cui le due popolazioni

tendono al punto di equilibrio situato in (x,y) = (𝑎𝜀2+𝑏𝑐

𝑏𝑑+𝜀2𝜀1,

𝑎𝑑−𝑐𝜀1

𝑏𝑑+𝜀2𝜀1).

Un’altro punto di equilibrio è il punto (x,y)=(0,0), ma si può facilmente vedere dal grafico sottostante che

esso non è di equilibrio stabile, infatti si può osservre che anche con un condizioni iniziali vicine al punto

(0,0), (in questo caso è posto x0=y0=0.0001) la traiettoria non converge a questo punto ma al precedente.

Questo conferma l’instabilità del punto (0,0).

Lo stesso ragionamento può essere adottato anche per il quarto punto di equilibrio, infatti si può vedere dal

seuente grafico che anche partendo molto vicino al punto i equilibrio (400,0) con condizioni iniziali

(400.0001 , 0.0001), esso non converge al punto di equilbrio in questione, ma al primo punto studiato:

Il terzo punto di equilibrio invece non ha senso studarlo perché non ha significato biologico (poiché richiede

un numero di predatori negativo).

Di seguito è riportato il grafico delle popolazioni di prede e predatori in funzione del tempo: possiamo

vedere che entrambe le popolazioni oscillano, con oscillazioni sempre più smorzate. L’apiezza di queste

oscillazioni non rimane costante ma tende a diminuire , come si può evincere dal grafico entrambe le

popolazioni tendo al

In quest’ultimo caso particolare a = 4, b = 2, c = 3, d = 1;

I punti di equilibrio sono (con 𝜀2 = 0.0001 e 𝜀1 = 0.0001) :

(𝑥, 𝑦) = (3.0002, 1.99985)

(𝑥, 𝑦) = (0,0);

(𝑥, 𝑦) = (0, − 30000 );

(𝑥, 𝑦) = (40000,0).

Invece sono (con 𝜀2 = 0.01 e 𝜀1 = 0.01) :

(𝑥, 𝑦) = (3.019849, 1.984901)

(𝑥, 𝑦) = (0,0);

(𝑥, 𝑦) = (0, − 300);

(𝑥, 𝑦) = (400,0).

È interessante anche vedere come si comporta il sistema per valori di 𝜀2 e 𝜀1 più alti (𝜀2 = 1 e 𝜀1 = 2),

dello stesso ordine di grandezza degli altri coefficienti a,b,c,d .

In questo caso la competizione tra le volpi è così forte a tal punto da estinguere le volpi stesse, mentre quella

tra le prede è tale da mantenere la popolazione costante nel tempo una volta arrivata al punto di equilibrio.

Infatti all’inizio sia le prede che i predatori diminuiscono, ma, mentre quella dei predatori diminuisce fino ad

estinguersi, ad un certo punto il numero delle prede inverte carattere ed inizia a risalire arrivando al punto di

equilibrio (in questo caso (3,0)). Con questa configurazione di parametri il punto di equilibrio stabile è

(𝑥, 𝑦) = (𝑎

𝜀1, 0) e non più (x,y) = (

𝑎𝜀2+𝑏𝑐

𝑏𝑑+𝜀2𝜀1,

𝑎𝑑−𝑐𝜀1

𝑏𝑑+𝜀2𝜀1), poiché in questo caso quest’ultimo punto ha coordinata

y negativa (quindi non ha senso biologico). Questo dimostra che nel caso precedente il punto (𝑥, 𝑦) =

(𝑎

𝜀1, 0) era instabile, mentre in questo caso risulta stabile, quindi possiamo affermare che esso è i realtà un

punto di sella.

Di seguito il grafico dell’andamento delle popolazioni in funzione del tempo.

In quest’ultimo caso particolare a = 6, b = 3, c = 4, d = 1, 𝜀2 = 1 e 𝜀1 = 2;

I punti di equilibrio sono:

(𝑥, 𝑦) = (3.6, −0.4)

(𝑥, 𝑦) = (0,0);

(𝑥, 𝑦) = (0, − 4 );

(𝑥, 𝑦) = (3,0).

MODELLO SUPER-PREDATORE (o con pesca)

In natura la distinzione tra prede e predatori non è sempre netta, i predatori, ad esempio, possono essere a

loro volta cacciati da una specie più forte. Spesso è utile tenere conto della presenza di un predatore esterno

in grado di cacciare entrambe le popolazioni dell'ecosistema.

Una situazione concreta è quella della pesca; durante questa pratica l'uomo caccia indistintamente due specie

diverse di pesci: una predatrice e l'altra predata

Per modellizare questa situazione è stato studiato questo sistema:

{

𝑑𝑥

𝑑𝑡= 𝑥(𝑎 − 𝑏𝑦 − 𝑒)

𝑑𝑦

𝑑𝑡= 𝑦(−𝑐 + 𝑑𝑥 − 𝑒)

Infatti possiamo immaginare che la pesca comporti il prelievo di individui in modo additivo e proporzionale

ad un parametro costante (𝒆 𝒐 𝒇), necessario a descrivere la pressione dell'uomo sul sistema.

I punti di equilibrio sono:

(𝑥, 𝑦) = (0, 0)

(𝑥, 𝑦) = (𝑐+𝑓

𝑑,

𝑎−𝑒

𝑏) ;

La stabilità dei punti di equilibrio può essere determinata linearizzando il sistema utilizzando la matrice

jacobiana del modello preda-predatore:

𝐽(𝑥, 𝑦) = [𝑎 − 𝑏𝑦 − 𝑒 −𝑏𝑥

𝑑𝑦 𝑑𝑥 − 𝑐 − 𝑒]

La matrice Jacobiana calcolata nel primo punto di equilibrio (0,0) è:

𝐽(0,0) = [𝑎 − 𝑒 0

0 −𝑐 − 𝑓]

Gli autovalori di questa matrice sono 𝜆1 = 1 = 𝜆2. Dato che almeno un autovalore della matrice J ha parte

reale positiva, allora il punto (0,0) è instabile.

La matrice jacobiana calcolata nel secondo punto di equilibrio è:

𝐽(𝑐 + 𝑒

𝑑,𝑎 − 𝑒

𝑏) = [

0 −𝑏(𝑐 + 𝑓)

𝑑𝑑(𝑎 − 𝑒)

𝑏0

]

Gli autovalori di questa matrice sono:

𝜆1 = 𝑖√(𝑐 + 𝑓)(𝑎 − 𝑒) e 𝜆1 = −𝑖√(𝑐 + 𝑓)(𝑎 − 𝑒)

Ci ritoviamo nello stesso caso del sistema lotka volterra classico in cui la parte reale degli autovalori è nulla

e perciò siamo in un caso particolare in cui non si possono trarre conclusioni dall'analisi lineare dela matrice

jacobiana. Come nel caso precedente il sistema ammette una costante del moto H(x,y) e le curve di livello di

questa funzione sono traiettorie chiuse che circondano il punto fisso, quindi, per lo stesso identico

ragionamento (vedi modello classico lotka volterrana sopra descirtto), possiamo dire si tratti di un punto di

equilibrio stabile.

Confrontando come in figura i ritratti in fase dei due sistemi (Lotka-Volterra classio e con pesca) e in

particolare le coordinate dei punti di equilibrio stabile di ciascun sistema si osserva subito che una pesca

moderata ha effetti che favoriscono la popolazione delle prede e sfavoriscono quella dei predatori: infatti

(𝑥, 𝑦) = (𝑐+𝑓

𝑑,

𝑎−𝑒

𝑏) 𝑝𝑢𝑛𝑡𝑜 𝑑𝑖 𝑒𝑞𝑢𝑖𝑙𝑖𝑏𝑟𝑖𝑜 𝑠𝑡𝑎𝑏𝑖𝑙𝑒 𝑐𝑜𝑛 𝑝𝑒𝑠𝑐𝑎;

(𝑥 , 𝑦 ) = (𝑐

𝑑,

𝑎

𝑏) punto di equilibrio stabile modello classico.

Per mostrare questo effetto si sono scelti dei parametri fissi (a = 2, b = 1, c = 2, d = 1) e si sono fatti variare i

parametri e-f relativi alla pesca.

Nel configurazione 1 del grafico sotto riportato sono stati scelti e=0.4 , f=0.6;

Nel configurazione 2 del sono stati scelti e=1.2 , f=1.8;

Nel configurazione 3 del sono stati scelti e=1.6 , f=2.4;

I punti di equilibrio delle diverse configurazioni sono:

P1 = (1.6, 1.6);

P2 = (2.8, 0.8);

P3= (3.4, 0.4).

Da ciò si può vedere che il punto di equilibrio stabile aumenta nella coordinata y e diminuisce nella x,

aumentando così all’equilibrio il numero di prede e diminuendo quelle di predatori.

È stata poi testata una configurazione limite impostando e=2 e f=3. In questo caso le prede non si

estinguono a differenza dei predatori. Osserviamo questo comportamento nel fictostante: si osserva la

popolazione di prede che, dopo un transitorio calo iniziale, si stabilizza su un valore costante, mentre i

predatori si estinguono esponenzialmente. (a,b,c,d sono rimasti invariati)

Oltrepassata questa soglia si ha l'estinzione di entrambe le specie in maniera esponenziale. All'aumentare dei

parametri e-f le due specie tendono all'estinzione più rapidamente.

Grafico che rappresenta i predatori in funzione delle prede: sono rappresentate 3 popolazioni in 3 colori

diversi; andando dal rosso al blu i coefficienti e-f aumentano.

Grafico delle popolazioni di prede in funzione del tempo: si veda come all’aumentare dei coefficienti e-f le

popolazioni si estinguono molto più rapidamente.

Come nel grafico precedente, anche questo grafico, che rappresenta le popolazioni dei predatori in funzione

del tempo all’aumentare dei coefficienti e-f le popolazioni si estinguono molto più rapidamente.

VERIFICA ALGORITMO DI INTEGRAZIONE

È qui rappresentato l’andamento della costante H del moto per una determinata combinazione di coefficienti

e condizioni iniziali (a=4, b=2, c=3, d=1, dt=0.0001, x0=8, y0=5), che si mantiene costante a meno di

piccolissime oscillazioni irregolari di ampiezza media dell’ordine di 10-13.

Come si può vedere dal grafico successivo, anche l’errore su H definito come 𝐻(𝑡)−𝐻0

𝐻0 presenta piccolissime

oscillazioni di ampiezza media dell’ordine di 10-14.

Nell’ultimo grafico è studiato l’andamento dell’errore sulla costante del moto H per diversi dt. Tale grafico

serve per verificare l’ordine d’integrazione dell’algoritmo scelto. Per Rounge Kutta di ordine quattro

l’andamento aspettato è una parabola di ordine 4 poiché di quarto ordine. Per rendere lo studio più semplice

si sono graficati i logaritmi in base 10 delle grandezze (dt e 𝐻(𝑡)−𝐻0

𝐻0 ) così da poter avere un andamento

lineare, in particolare è necessario verificare che il coefficiente angolare della retta sia 4. Per studiare l’errore

di integrazione è stato preso il valore di H per ogni scelta del dt (il dt veniva ogni volta dimezzato) sempre

allo stesso tempo t=1. La retta ottenuta, tramite un fit lneare ha coefficiente angolare m=4.012±0.016 ed

intercetta q=3.26±0.10.

Il risultato ottenuto è soddisfacente, il coefficiente angolare è compatibile con quello aspettato.

Nel box sottostante sono riportati i dati riguardati lo studio di integrazione del metodo utilizzato, con i

parametri specificati all’inizio).

t dt (H(T)-H0)/H0)

1.0000000000 0.03125000000000 0.00002670462687322557

1.0000000000 0.01562500000000 0.00000149749755925359

1.0000000000 0.00781250000000 0.00000008829112851873

1.0000000000 0.00390625000000 0.00000000535308865626

1.0000000000 0.00195312500000 0.00000000032941731647

1.0000000000 0.00097656250000 0.00000000002042772901

1.0000000000 0.00048828125000 0.00000000000127172823

1.0000000000 0.00024414062500 0.00000000000008291341

1.0000000000 0.00012207031250 0.00000000000000600580

Di seguito è riportato il file con i dati del fit lineare.

*******************************************************************************

Tue Dec 30 12:28:20 2014

FIT: data read from "LVerr.dat" u (log($2)):(log($5))

format = x:z

#datapoints = 9

residuals are weighted equally (unit weight)

function used for fitting: f(x)

fitted parameters initialized with current variable values

Iteration 0

WSSR : 2721.44 delta(WSSR)/WSSR : 0

delta(WSSR) : 0 limit for stopping : 1e-05

lambda : 4.64326

initial set of free parameter values

a = 1

b = 1

After 5 iterations the fit converged.

final sum of squares of residuals : 0.0513647

rel. change during last iteration : -2.90445e-14

degrees of freedom (FIT_NDF) : 7

rms of residuals (FIT_STDFIT) = sqrt(WSSR/ndf) : 0.0856611

variance of residuals (reduced chisquare) = WSSR/ndf : 0.00733782

Final set of parameters Asymptotic Standard Error

======================= ==========================

a = 3.2574 +/- 0.1035 (3.179%)

b = 4.01199 +/- 0.01595 (0.3977%)

correlation matrix of the fit parameters:

a b

a 1.000

b 0.961 1.000

LISTATO

#include<stdio.h> // Equazione Lortka-Volterra RK4

#include<stdlib.h>

#include<math.h>

#include<time.h>

void LotkaVolterraStandard(double a,double b,double c, double d,double dt,double x0,double y0);

void LotkaVolterraSpecial(double a, double b, double c,double d, double e1, double e2,double

dt,double x0,double y0);

void LotkaVolterraPesca(double a,double b,double c, double d,double e,double f,double dt,double

x0,double y0);

void Errore(double a,double b,double c, double d,double dt,double x0,double y0);

void FunzioneH(double a,double b,double c, double d,double dt,double x0,double y0);

void Inizializzazione();

int main(){

Inizializzazione();

}

void LotkaVolterraStandard(double a,double b,double c, double d,double dt,double x0,double y0){

double k1,k2,k3,k4,l1,l2,l3,l4,X,Y,t,H;

int i=0;

double J0[1][1],J1[1][1];

J0[0][0]=a;

J0[0][1]=0;

J0[1][0]=0;

J0[1][1]=-c;

J1[0][0]=0;

J1[0][1]=-(b*c)/d;

J1[1][0]=(a*d)/b;

J1[1][1]=0;

printf("I punti di equilibrio sono:\n1)(x,y)=(0,0)\n2)(x,y)=(%lf,%lf)\n\n",c/d,a/b);

printf("La matrice jacobiana calcolata nel primo punto di equilibrio (0,0) è: \n\n %lf

%lf \n",J0[0][0],J0[0][1]);

printf("J(0,0)=\n");

printf(" %lf %lf \n\n",J0[1][0],J0[1][1]);

printf("La matrice jacobiana calcolata nel secondo punto di equilibrio (%.2lf,%.2lf) è: \n\n

%lf %lf \n",c/d,a/b,J1[0][0],J1[0][1]);

printf("J(%.2lf,%.2lf)=\n",c/d,a/b);

printf(" %lf %lf \n",J1[1][0],J1[1][1]);

FILE *fp;

fp=fopen("LVstd.dat","w+");

fprintf(fp,"#x0 y0 t

H\n");

t=0;

for(i=0;i<=1000000;i++){

H=-(c*log(x0)-d*x0-b*y0+a*log(y0));

k1=x0*(a-b*y0);

l1=y0*(-c+d*x0);

k2=(x0+0.5*dt*k1)*(a-b*(y0+0.5*dt*l1));

l2=(y0+0.5*dt*l1)*(-c+d*(x0+0.5*dt*k1));

k3=(x0+0.5*dt*k2)*(a-b*(y0+0.5*dt*l2));

l3=(y0+0.5*dt*l2)*(-c+d*(x0+0.5*dt*k2));

k4=(x0+dt*k3)*(a-b*(y0+dt*l3));

l4=(y0+dt*l3)*(-c+d*(x0+dt*k3));

X=x0+(1./6)*dt*(k1+2*k2+2*k3+k4);

Y=y0+(1./6)*dt*(l1+2*l2+2*l3+l4);

fprintf(fp,"%.10lf %.10lf %lf

%.20lf\n",x0,y0,t,H);

x0=X;

y0=Y;

t=t+dt;

}

fclose(fp);

}

void LotkaVolterraSpecial(double a, double b, double c,double d, double e1, double e2,double

dt,double x0,double y0){

double m1,m2,m3,m4,n1,n2,n3,n4,X2,Y2,t2,H,X,Y;

int h=0;

double J0[1][1],J1[1][1];

X=(a*e2+b*c)/(b*d+e1*e2);

Y=(a*d-c*e1)/(b*d+e1*e2);

J0[0][0]=a;

J0[0][1]=0;

J0[1][0]=0;

J0[1][1]=-c;

J1[0][0]=-e1*X;

J1[0][1]=-b*X;

J1[1][0]=d*Y;

J1[1][1]=-e2*Y;

printf("I punti di equilibrio sono:\n1)(x,y)=(0,0)\n2)(x,y)=(%lf,%lf)\n3) (x,y)=(%lf,%lf)\n4)

(x,y)=(%lf,%lf)\n\n",(a*e2+b*c)/(b*d+e1*e2),(a*d-c*e1)/(b*d+e1*e2),0.,-c/e2,a/e1,0.);

printf("La matrice jacobiana calcolata nel primo punto di equilibrio (0,0) è: \n\n

%lf %lf \n",J0[0][0],J0[0][1]);

printf("J(0,0)=\n");

printf(" %lf %lf \n\n",J0[1][0],J0[1][1]);

printf("La matrice jacobiana calcolata nel secondo punto di equilibrio (%lf,%lf) è: \n\n

%lf %lf \n",X,Y,J1[0][0],J1[0][1]);

printf("J(%lf,%lf)=\n",X,Y);

printf(" %lf %lf \n",J1[1][0],J1[1][1]);

FILE *fq;

fq=fopen("LVsp.dat","w+");

t2=0;

fprintf(fq,"#x0 y0 t

H\n");

for(h=0;h<=100000;h++){

H=c*log(x0)-d*x0-b*y0+a*log(y0);

m1=x0*(a-b*y0-e1*x0);

n1=y0*(-c+d*x0-e2*y0);

m2=(x0+0.5*dt*m1)*(a-b*(y0+0.5*dt*n1)-e1*(x0+0.5*dt*m1));

n2=(y0+0.5*dt*n1)*(-c+d*(x0+0.5*dt*m1)-e2*(y0+0.5*dt*n1));

//printf("%lf %lf\n",m2,n2);

m3=(x0+0.5*dt*m2)*(a-b*(y0+0.5*dt*n2)-e1*(x0+0.5*dt*m2));

n3=(y0+0.5*dt*n2)*(-c+d*(x0+0.5*dt*m2)-e2*(y0+0.5*dt*n2));

m4=(x0+dt*m3)*(a-b*(y0+dt*n3)-e1*(x0+dt*m3));

n4=(y0+dt*n3)*(-c+d*(x0+dt*m3)-e2*(y0+dt*n3));

X2=x0+(1./6)*dt*(m1+2*m2+2*m3+m4);

Y2=y0+(1./6)*dt*(n1+2*n2+2*n3+n4);

fprintf(fq,"%lf %lf %lf %lf\n",x0,y0,t2,H);

x0=X2;

y0=Y2;

t2=t2+dt;

}

fclose(fq);

}

void LotkaVolterraPesca(double a,double b,double c, double d,double e,double f,double dt,double

x0,double y0){

double m1,m2,m3,m4,n1,n2,n3,n4,X2,Y2,t2,H,X,Y;

int h=0;

double J0[1][1],J1[1][1];

J0[0][0]=a-e;

J0[0][1]=0;

J0[1][0]=0;

J0[1][1]=-c-f;

J1[0][0]=0;

J1[0][1]=(-b*(c+f))/d;

J1[1][0]=(d*(a-e))/b;

J1[1][1]=0;

printf("I punti di equilibrio sono:\n1)(x,y)=(0,0)\n2)(x,y)=(%lf,%lf)\n\n",(c+f)/d,(a-e)/b);

printf("La matrice jacobiana calcolata nel primo punto di equilibrio (0,0) è: \n\n

%lf %lf \n",J0[0][0],J0[0][1]);

printf("J(0,0)=\n");

printf(" %lf %lf \n\n",J0[1][0],J0[1][1]);

printf("La matrice jacobiana calcolata nel secondo punto di equilibrio (%lf,%lf) è: \n\n

%lf %lf \n",(c+f)/d,(a-e)/b,J1[0][0],J1[0][1]);

printf("J(%lf,%lf)=\n",(c+f)/d,(a-e)/b);

printf(" %lf %lf \n",J1[1][0],J1[1][1]);

FILE *fu;

fu=fopen("LVfish.dat","w+");

t2=0;

fprintf(fu,"#x0 y0 t

H\n");

for(h=0;h<=100000;h++){

H=-((c+f)*log(x0)-d*x0-b*y0+(a-e)*log(y0));

m1=x0*(a-b*y0-e);

n1=y0*(-c+d*x0-f);

m2=(x0+0.5*dt*m1)*(a-b*(y0+0.5*dt*n1)-e);

n2=(y0+0.5*dt*n1)*(-c+d*(x0+0.5*dt*m1)-f);

//printf("%lf %lf\n",m2,n2);

m3=(x0+0.5*dt*m2)*(a-b*(y0+0.5*dt*n2)-e);

n3=(y0+0.5*dt*n2)*(-c+d*(x0+0.5*dt*m2)-f);

m4=(x0+dt*m3)*(a-b*(y0+dt*n3)-e);

n4=(y0+dt*n3)*(-c+d*(x0+dt*m3)-f);

X2=x0+(1./6)*dt*(m1+2*m2+2*m3+m4);

Y2=y0+(1./6)*dt*(n1+2*n2+2*n3+n4);

fprintf(fu,"%lf %lf %lf %lf\n",x0,y0,t2,H);

x0=X2;

y0=Y2;

t2=t2+dt;

}

fclose(fu);

}

void Errore(double a,double b,double c, double d,double dt,double x0,double y0){

FILE *fs;

fs=fopen("LVerr.dat","w+");

fprintf(fs,"#t dt H H0

|(H-H0)/H0|\n");

double k1,k2,k3,k4,l1,l2,l3,l4,X,Y,t,H,H0,xi,yi;

H0=-(c*log(x0)-d*x0-b*y0+a*log(y0));

//printf("%lf\n",log(2.7182818284));

xi=x0;

yi=y0;

do{

x0=xi;

y0=yi;

t=0;

while(t<100){

//printf("%lf\n",t);

H=-(c*log(x0)-d*x0-b*y0+a*log(y0));

k1=x0*(a-b*y0);

l1=y0*(-c+d*x0);

k2=(x0+0.5*dt*k1)*(a-b*(y0+0.5*dt*l1));

l2=(y0+0.5*dt*l1)*(-c+d*(x0+0.5*dt*k1));

k3=(x0+0.5*dt*k2)*(a-b*(y0+0.5*dt*l2));

l3=(y0+0.5*dt*l2)*(-c+d*(x0+0.5*dt*k2));

k4=(x0+dt*k3)*(a-b*(y0+dt*l3));

l4=(y0+dt*l3)*(-c+d*(x0+dt*k3));

X=x0+(1./6)*dt*(k1+2*k2+2*k3+k4);

Y=y0+(1./6)*dt*(l1+2*l2+2*l3+l4);

x0=X;

y0=Y;

t=t+dt;

if(t>2-0.5*dt && t<2+0.5*dt){

//printf(" %lf %.14lf %lf %lf

%.20lf\n",t,dt,H0,H0,fabs(((H-H0)/H0)));

fprintf(fs," %lf %.14lf %.10lf %.10lf

%.20lf\n",t,dt,H,H0,fabs(((H-H0)/H0)));

}

}

dt=dt/2;

}while(dt>0.0001);

fclose(fs);

}

void FunzioneH(double a,double b,double c, double d,double dt,double x0,double y0){

double k1,k2,k3,k4,l1,l2,l3,l4,X,Y,t,H,H0,k,l;

double tmp1,tmp2;

tmp1=x0;

tmp2=y0;

int i=0;

FILE *ft;

ft=fopen("LVfun.dat","w+");

fprintf(ft,"#x0 y0 H0\n");

for(k=1;k<=tmp1;k++){

printf("%lf\n",k);

x0=k;

for(l=1;l<=tmp2;l++){

H0=-(c*log(x0)-d*x0-b*l+a*log(l));

t=0;

y0=l;

for(i=0;i<10000;i++){

H=c*log(x0)-d*x0-b*y0+a*log(y0);

k1=x0*(a-b*y0);

l1=y0*(-c+d*x0);

k2=(x0+0.5*dt*k1)*(a-b*(y0+0.5*dt*l1));

l2=(y0+0.5*dt*l1)*(-c+d*(x0+0.5*dt*k1));

k3=(x0+0.5*dt*k2)*(a-b*(y0+0.5*dt*l2));

l3=(y0+0.5*dt*l2)*(-c+d*(x0+0.5*dt*k2));

k4=(x0+dt*k3)*(a-b*(y0+dt*l3));

l4=(y0+dt*l3)*(-c+d*(x0+dt*k3));

X=x0+(1./6)*dt*(k1+2*k2+2*k3+k4);

Y=y0+(1./6)*dt*(l1+2*l2+2*l3+l4);

fprintf(ft,"%lf %lf %lf\n",x0,y0,H0);

x0=X;

y0=Y;

t=t+dt;

}

}

}

fclose(ft);

}

void Inizializzazione(){

double dt,x0,y0,e1,e2;

double a,b,c,d,f1,f2,e,f;

int q,w,i;

double seed;

seed=time(0);

srand48(seed);

printf("\nQuale modello vuoi studiare?\n\n1)Modello classico...........................(premere

1)\n2)Modello con competizione intraspecie.......(premere 2)\n3)Modello con super-

predatore................(premere 3)\n\n");

scanf("%d",&q);

printf("\n");

printf("In che modo vuoi procedere?\n\n1)Inserimento numeri da default.......(premere

1)\n2)Inseriemento numeri da tastiera.....(premere 2)\n3)Inserimento numeri

random...........(premere 3)\n\n");

scanf("%d",&w);

printf("\n");

if(w==1){

dt=0.001;

x0=8;

y0=5;

if(q==1){

a=4;

b=2;

c=3;

d=1;

LotkaVolterraStandard(a,b,c,d,dt,x0,y0);

FunzioneH(a,b,c,d,0.001,x0,y0);

}

if(q==2){

a=6;

b=3;

c=4;

d=1;

f1=0.01;

f2=0.01;

LotkaVolterraSpecial(a,b,c,d,f1,f2,dt,x0,y0);

}

if(q==3){

a=4;

b=2;

c=3;

d=1;

e=2;

f=1;

LotkaVolterraPesca(a,b,c,d,e,f,dt,x0,y0);

FunzioneH(a-e,b,c+f,d,0.001,x0,y0);

}

Errore(a,b,c,d,0.03,x0,y0);

}

if(w==2){

printf("dimmi il passo di integrazione (dt)\n");

scanf("%lf",&dt);

printf("dimmi il numero iniziale di prede (x0)\n");

scanf("%lf",&x0);

printf("dimmi il numero iniziale di predatori (y0)\n");

scanf("%lf",&y0);

printf("dimmi a\n");

scanf("%lf",&a);

printf("dimmi b\n");

scanf("%lf",&b);

printf("dimmi c\n");

scanf("%lf",&c);

printf("dimmi d\n");

scanf("%lf",&d);

if(q==1){

LotkaVolterraStandard(a,b,c,d,dt,x0,y0);

FunzioneH(a,b,c,d,dt,x0,y0);

}

if(q==2){

printf("dimmi f1\n");

scanf("%lf",&f1);

printf("dimmi f2\n");

scanf("%lf",&f2);

LotkaVolterraSpecial(a,b,c,d,f1,f2,dt,x0,y0);

}

if(q==3){

printf("dimmi e\n");

scanf("%lf",&e);

printf("dimmi f\n");

scanf("%lf",&f);

LotkaVolterraPesca(a,b,c,d,e,f,dt,x0,y0);

FunzioneH((a-e),b,(c+f),d,0.001,x0,y0);

}

Errore(a,b,c,d,0.03,x0,y0);

}

if(w==3){

a=20*drand48();

b=20*drand48();

c=20*drand48();

d=20*drand48();

x0=(int)(20*drand48());

y0=(int)(20*drand48());

dt=0.001;

printf("I numeri estratti sono:\n\na=%lf b=%lf c=%lf d=%lf x0=%lf

y0=%lf\n\n",a,b,c,d,x0,y0);

if(q==1){

LotkaVolterraStandard(a,b,c,d,dt,x0,y0);

FunzioneH(a,b,c,d,dt,x0,y0);

}

if(q==2){

f1=20*drand48();

f2=20*drand48();

printf("f1=%lf f2=%lf \n\n",f1,f2);

LotkaVolterraSpecial(a,b,c,d,f1,f2,dt,x0,y0);

}

if(q==3){

e=20*drand48();

f=20*drand48();

LotkaVolterraPesca(a,b,c,d,e,f,dt,x0,y0);

FunzioneH(a-e,b,c+f,d,0.001,x0,y0);

}

Errore(a,b,c,d,0.03,x0,y0);

}

}

//FINE Equazione Lotka-Volterra RK4