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
Top Related