Corso di Sistemi caotici per la protezione dell’informazione...

Università degli Studi di Genova

Dipartimento di Ingegneria Biofisica ed Elettronica

Corso di Sistemi caotici per la protezione

dell’informazione 1

Docente: Prof. S. Cincotti

Studio del modello di Rosenzweig-MacArthur con relativo sviluppo di un'interfaccia

di simulazione in ambiente Matlab

Manuale d’uso

Enrico Costa Paolo Spiga

Genova, 18 maggio 2004

Enrico Costa, Paolo Spiga Studio del modello di Rosenzweig-MacArthur con relativo sviluppo di un'interfaccia di simulazione in ambiente Matlab

2

INDICE

1.0 INTRODUZIONE pag. 3

1.1 Modello di Rosenzweig-MacArthur pag. 5

1.1.1 Analisi dei punti fissi pag. 6

1.2 Considerazioni generali ed approcci allo studio di biforcazione del modello in letteratura

pag. 8 2.0 INTERFACCIA DI SIMULAZIONE

pag. 11 2.1 Descrizione generale dell'interfaccia grafica

pag. 12 2.2 Simulazione dell'evoluzione temporale del modello

di Rosenzweig-MacArthur: verifica delle tipologie di comportamento del sistema dinamico

pag. 13 2.2.1 Estinzione del superpredatore pag. 13 2.2.2 Coesistenza stazionaria pag. 17 2.2.3 Coesistenza ciclica in bassa frequenza pag. 20 2.2.4 Coesistenza caotica pag. 22 2.2.5 Coesistenza ciclica in alta frequenza pag. 24

3.0 ANALISI DI BIFORCAZIONE pag. 26

3.1 Generalità sull’utilizzo di Content pag. 27

3.2 Risultati dell’analisi pag. 28

4.0 CONCLUSIONI pag. 36

Riferimenti bibliografici pag. 37


3

1.0 INTRODUZIONE Lo studio delle dinamiche di popolazione è oggi di grande interesse in diversi ambiti scientifici; numerosi e differenti sono, infatti, gli obbiettivi che si possono raggiungere avendo a disposizione un modello matematico che descriva con un buon grado di accuratezza i sistemi dinamici presenti in natura e si possono pertanto individuare molteplici approcci verso questi temi nella letteratura di discipline in generale non strettamente connesse tra loro. Anche se l’ingegneria non è direttamente interessata al comportamento dei sistemi biologici, tuttavia non si può trascurare che i comportamenti caotici, riscontrabili in natura, destano sempre maggiore interesse nell’ambito della realizzazione di strumenti elettronici ed informatici per la vita quotidiana di tutti, non solo degli addetti ai lavori: si pensi, per citare alcuni esempi tutt’altro che esaustivi dell’argomento, al problema di generazione di chiavi e di numeri casuali, oppure alla sincronizzazione di sistemi caotici per la trasmissione sicura dell’informazione, nell’ambito della crittografia; insomma, lo studio dei sistemi dinamici è oggetto di studi trasversali rispetto alle varie discipline che portano il proprio peculiare contributo ad una conoscenza più precisa dei fenomeni di tipo caotico. Le dinamiche di popolazione studiano l’andamento nel tempo della densità degli esemplari di una o più specie animali o vegetali, considerando l’ambiente in cui si sviluppano e l’interazione con altre specie o fenomeni esterni; in particolare, si cercano modelli matematici, più o meno complessi in funzione del numero di specie osservate e di fenomeni considerati determinanti per l’evoluzione del sistema dinamico, che consentano di fare previsioni sui possibili interventi esterni per prevedere e guidare il comportamento delle popolazioni coinvolte verso situazioni ottimali o di interesse. A questo proposito si può notare che i comportamenti dei sistemi naturali evolvono spesso in maniera assai complessa; in particolare le dinamiche irregolari osservabili in molte catene alimentari possono essere quelle di strani attrattori [2], che talvolta sono compresenti e giustificano così, in condizioni particolari, cambiamenti repentini e radicali nel comportamento qualitativo del sistema [3], senza minare tuttavia la persistenza delle specie interessate nell’ecosistema considerato: il principio generale “optimality implies chaos”, richiamato in [1,2] e presente in numerose aree della biologia, dimostra così la sua validità anche in questo caso. Inoltre è necessario sottolineare che la presenza di regimi ciclici vicini a comportamenti caotici non garantisce che un arricchimento sconsiderato della catena alimentare, attraverso l’aumento della risorsa primaria del livello trofico più basso, implichi automaticamente la persistenza di tutta la catena alimentare, in accordo con le argomentazioni addotte da Rosenzweig,


4

per catene alimentari ditrofiche, circa il “paradosso dell’arricchimento” [4] che mettono in guardia dall’intraprendere semplicistiche iniziative per il controllo degli ecosistemi. Dobbiamo notare che, per descrivere realisticamente ed in maniera sufficientemente precisa un insieme di popolazioni che interagiscono in un ambiente, non è possibile limitarsi all’utilizzo di funzioni lineari proprio a causa dell’intrinseca complessità, che talvolta diventa “caoticità”, dei sistemi naturali [2]; tuttavia il necessario inserimento di non-linearità può determinare un tale aumento della complessità del modello matematico da rendere necessario l’ausilio di software, più o meno specifico per l’analisi dei sistemi dinamici, che fornisca dati più ampi e dettagliati di quelli ottenibili mediante lo studio analitico classico, ferma restando la necessità di un’interpretazione critica dei risultati per la definizione del comportamento complessivo del sistema considerato. Il modello matematico utilizzato per la descrizione del sistema dinamico di interesse può essere descritto in tempo continuo, attraverso un sistema di equazioni differenziali, oppure in tempo discreto, mediante mappe iterate, ed in entrambi i casi l’ordine del sistema corrisponde al numero di specie che vengono studiate; nel primo caso si definisce la derivata della concentrazione di ciascuna specie in funzione delle densità di tutte le popolazioni costituenti il sistema, mentre nella seconda rappresentazione si fornisce l’espressione delle concentrazioni delle diverse specie rispetto ai loro valori al passo precedente. E’ necessario ricordare che, a differenza della rappresentazione in tempo discreto, l’utilizzo di modelli in tempo continuo prevede la possibilità di individuare comportamenti caotici solamente per sistemi di ordine tre o superiore; per questo motivo il nostro studio ha preso in considerazione il modello tritrofico di Rosenzweig-MacArthur [1, 2, 3] che descrive una catena alimentare costituita da una preda, un predatore ed un superpredatore interagenti tra loro ed inseriti nello stesso ambiente.


5

1.1 Modello di Rosenzweig - MacArthur Come già indicato nell’introduzione, il modello di Rosenzweig-MacArthur si propone di studiare il comportamento di una catena alimentare tritrofica costituita da una preda (o risorsa primaria), un predatore ed un superpredatore. In particolare questo modello a tempo continuo prevede la presenza di una preda con comportamento logistico e di un predatore ed un superpredatore di tipo Holling II [1, 2] descritti dal sistema di equazioni

1 1 2 21

2 1

3 32 2 12 2 2

2 1 3 2

3 3 23 3 3

3 2

1dx x a x

x rdt K b x

a xdx a xx e d

dt b x b x

dx a xx e d

dt b x

= − − + = − − + + = − +

(1)

dove 1x , 2x e 3x sono le densità rispettivamente di preda, predatore e superpredatore, r e K sono la natalità e la capacità di autosostentamento della preda, mentre ia , ib , ie e id sono rispettivamente la velocità massima di predazione, la costante di saturazione parziale, l’efficienza e la mortalità del predatore ( 2i = ) e del superpredatore ( 3i = ). Analizzando semplicemente le equazioni del modello e la loro tipologia, dobbiamo innanzitutto evidenziare che, per come è stato definito il modello, tutti i parametri si devono intendere positivi ( K , 2b e 3b devono essere strettamente positivi; gli altri parametri possono annullarsi, causando tuttavia la degenerazione del sistema dovuta a fenomeni di interazione trascurati) e possiamo quindi notare che: - la densità della preda tende in generale a crescere per valori inferiori alla sua

capacità di autosostentamento ed a decrescere con la densità del predatore; - la densità del predatore tende in generale a crescere con se stessa e con la

densità della preda, ma viene limitata e contrastata dalla sua mortalità e dalla densità del superpredatore;

- la densità del superpredatore tende in generale a crescere con se stessa e con la densità del predatore, ma viene limitata dalla sua mortalità.

Inoltre è bene sottolineare che non tutte le configurazioni dei parametri sono in generale accettabili come descrittive di catene alimentari realistiche; in particolare vogliamo evidenziare qui la condizione 222 dea > ( 333 dea > ) che garantisce che il predatore (il superpredatore) sopravviva per concentrazioni infinite della preda (del predatore), cioè che la concentrazione del consumatore


6

non si annulli proprio in corrispondenza dell’infinita disponibilità della sua risorsa [2, 3]. 1.1.1 Analisi dei punti fissi Ci interessiamo ora alla ricerca dei punti fissi (punti di equilibrio) [3], cioè tutti quei punti che soddisfano, per il sistema di equazioni (1) con tutti i parametri strettamente positivi, la condizione 0=dtdx i ; riportiamo tutte le soluzioni analitiche del sistema, indicando tuttavia esplicitamente le condizioni sul valore dei parametri per le quali esse appartengono all’ottante positivo dello spazio tridimensionale ( )321 ,, xxx e sono quindi significative per lo studio del comportamento dinamico della catena alimentare tritrofica considerata: q ( ) ( )0,0,0,,

*

1321 =xxx , ossia estinzione delle tre specie;

q ( ) ( )0,0,,,*

2321 Kxxx = , ossia estinzione di predatore e superpredatore e sopravvivenza della preda che attesta la sua densità sul valore della capacità di autosostentamento;

q ( ) ( )

−= 0,,,,

2

2*

3321 HKHKdre

Hxxx con 222

22

deadb

H−

= , ossia estinzione del

superpredatore e convivenza di preda e predatore sotto la condizione generale 222 dea > già descritta e quella particolare KH < (per 0=H e

Kdea

dbH =

−=

222

22 si ottengono rispettivamente ( )*

1321 ,, xxx e ( )*

2321 ,, xxx );

q ( ) ( )ZYXxxx ,,,,*

4,4321 =−+ con ( ) ( )

2

4 2222 Y

rKa

bKbKX

−+±−= ,

333

33

deadb

Y−

=

e

−

+= 2

2

22

3

3

1d

Xbea

Yd

eZ , ossia coesistenza delle tre specie sotto la

condizione generale 333 dea > già descritta e quelle particolari 222

22

deadb

X−

>

e ( ) 333

332

2

24

deadb

bK

Kar

−+> (per

333

33

2

2

deadb

ba

r−

> si ha una sola soluzione

“positiva” di X , mentre per 333

33

2

2

deadb

ba

r−

< e 2bK > le soluzioni sono due);

q esiste un ulteriore punto fisso, di coordinate

−

−− 333

332

333

33 ,,0dea

ebddea

db, privo

di significato fisico per il sistema dinamico da noi considerato, dal momento che in generale non sono accettabili valori negativi per variabili di stato che rappresentino la concentrazione di popolazioni costituenti una catena alimentare.


7

Per quanto riguarda la stabilità dei punti fissi, occorre calcolare gli autovalori della matrice jacobiana del sistema, che nella sua forma generale è

( )

( ) ( )

( )

−++

−+

−++

+−

+−−

=

323

2332

23

3333

2223

333

12

1222

12

2222

12

122

12

2221

0

0

02

dxb

xea

xb

xeba

dxb

xbaxbxea

xb

xebaxb

xa

xb

xbaKrx

r

J .

Se non è agevole la valutazione di J per ( )*

3321 ,, xxx e ( )*

4,4321 ,, −+xxx , tuttavia è

facile ricavare gli autovalori per ( )*

1321 ,, xxx e ( )*

2321 ,, xxx che valgono rispettivamente

0,,0

0,0,0

3232

2

222221

31321211

<−=−+

=<−=

<−=<−=>=

ddKbKea

r

ddr

λλλ

λλλ

;

in particolare si può notare che l’origine è globalmente instabile a causa di un

autovalore positivo, mentre l’equilibrio ( )0,0,K è stabile per 22

22 dKbKea

<+

.


8

1.2 Considerazioni generali ed approcci allo studio di biforcazione del modello in letteratura

Come si può evincere dai risultati dell’analisi per la ricerca dei punti fissi, le soluzioni trovate non permettono di individuare, per una qualsiasi configurazione dei parametri, una situazione di possibile convivenza stabile per le tre specie analizzate e pertanto non sono sistematicamente utilizzabili ai fini di un controllo esterno per la tutela di un sistema naturale riconducibile al modello di Rosenzweig-MacArthur; perciò, l’interesse dimostrato dagli studiosi [1, 2 e cit.] per questo modello si è orientato verso l’individuazione dei possibili interventi umani non tanto per garantire che le tre specie possano convivere in un equilibrio stabile, quanto per evitare che si estinguano gli esemplari anche di uno solo dei gradini della scala trofica. Inoltre bisogna sempre ricordare che qualunque tentativo per la definizione di modelli relativi ai sistemi naturali inevitabilmente trascura e minimizza i contributi e gli effetti di alcuni fattori esterni e pertanto è necessario considerare la possibilità di piccoli disallineamenti tra modello matematico e sistema reale: risulterebbe addirittura rischioso intraprendere azioni di intervento senza considerare questa possibilità ed è pertanto opportuno cercare di individuare nei diagrammi di biforcazione regioni dei parametri che presentino tipologie di comportamento che garantiscano la sopravvivenza di tutte le specie della catena alimentare. Ovviamente, per la definizione stessa del modello e come peraltro ci si può aspettare intuitivamente, la specie più a rischio di estinzione è quella del superpredatore e pertanto la possibilità di garantire la sua sopravvivenza fornisce una possibile indicazione sulle azioni da intraprendere per rafforzare la catena alimentare tritrofica nel suo complesso. In particolare gli studi da noi considerati [1, 2], che si sono posti come obbiettivo la massimizzazione delle densità delle tre specie, hanno individuato appunto in 3x l’indicatore di “buona salute” di tutta la catena ed hanno considerato le variazioni del comportamento qualitativo del sistema rispetto ai parametri r e K ; la scelta di questi due parametri, in particolar modo quella di K , è stata guidata dalla necessità di individuare la possibilità, successivamente confermata dallo studio di biforcazione, di controllare dall’esterno la catena alimentare mediante il nutrimento della preda, che influisce appunto sul valore di r e, in modo ancora più determinante, su quello di K [1]. Per orientarsi nello spazio dei parametri durante le simulazioni che abbiamo condotto sul modello di Rosenzweig-MacArthur è stato utilizzato il seguente diagramma tratto da [1, 2] che, attraverso alcune regioni individuate nel piano ( K , r ), indica diverse tipologie di comportamento per le tre specie della catena alimentare. In particolare il diagramma è stato ottenuto fissando gli altri parametri del sistema a valori frequentemente presenti in letteratura [1 e cit.]


9

( 352 =a , 312 =b , 12 =e , 1042 =d , 2013 =a , 213 =b , 13 =e e 10013 =d ) che sono stati utilizzati anche da noi in tutte le nostre simulazioni; segnaliamo tuttavia che esistono altre numerose e differenti configurazioni possibili per i parametri [3] che non rientrano, però, nella nostra trattazione.

In generale possiamo notare che per bassi valori di r il superpredatore si estingue a causa della ridotta rigenerazione della preda, mentre per valori di r superiori troviamo comportamenti più o meno complessi al variare del parametro K , che modella l’influenza del nutrimento della preda e quindi l’arricchimento di tutta la catena alimentare; in particolare, per un arricchimento continuo del ciclo alimentare (K crescente), la complessità del comportamento dinamico del sistema in primo luogo aumenta da stazionario a ciclico, passa poi alla coesistenza caotica e infine diminuisce, da caotico alla coesistenza ciclica: la coesistenza è dunque possibile, secondo diverse modalità, per intervalli non trascurabili dello spazio dei parametri. Non deve meravigliare la presenza di una regione, per quanto ridotta, di comportamento caotico, che la teoria prevede possibile nei sistemi di ordine tre o superiore. Come si può notare nel grafico, la regione di coesistenza caotica è piuttosto stretta ed allungata, delimitata a destra da una curva regolare ed a sinistra da una struttura frattale prodotta dalla biforcazione [1]; con un’analisi più approfondita si può scoprire la presenza di strani attrattori con forma a “tazzina da tè” (teacup), al confine con la regione di coesistenza ciclica a bassa frequenza, e a “tazzina da tè tagliata” (cut teacup), al confine con la regione di coesistenza ciclica ad alta frequenza, che manifestano la loro somiglianza con le forme tipiche degli attrattori delle adiacenti tipologie di comportamento.


10

Il precedente grafico è molto utile per individuare in maniera immediata il comportamento qualitativo del sistema dinamico al variare dei parametri r e K , ma è opportuno considerare il diagramma di biforcazione da cui esso deriva [1] per avere un punto di riferimento e di conferma per la nostra analisi sul modello.

Nello studio [1] condotto servendosi del programma AUTO [5], le linee continue si riferiscono a biforcazioni catastrofiche di attrattori, mentre quelle tratteggiate indicano biforcazioni non-catastrofiche di attrattori o biforcazioni di punti sella o di repulsori; in particolare i simboli riportati nel grafico hanno il seguente significato: - TE : transcritical of equilibria; - FE : fold of equilibria; - Hp : planar Hopf; - H� : supercritical Hopf; - H+ : subcritical Hopf; - TC : transcritical of cycles; - FC : fold of cycles. Nelle nostre considerazioni non ci soffermeremo in maniera puntuale su questi risultati, per i quali rimandiamo a [1], ma ci occuperemo di verificare ed ampliare l’analisi dei punti fissi illustrata nel precedente paragrafo.


11

2.0 INTERFACCIA DI SIMULAZIONE Come abbiamo già evidenziato nell’introduzione e nel primo capitolo, il modello di Rosenzweig-MacArthur presenta un consistente numero di parametri, ma nella nostra analisi considereremo solamente l’approccio, presente in letteratura [1, 2], di uno studio di biforcazione rispetto ai parametri r e K ; in particolare, in questo capitolo rivolgeremo la nostra attenzione alla verifica dei risultati di alcune simulazioni condotte in [6] circa le tipologie di comportamento del sistema dinamico mediante l’uso di un interfaccia grafica in ambiente MATLAB da noi sviluppata appositamente e di seguito brevemente descritta.


12

2.1 Descrizione generale dell’interfaccia grafica In questo paragrafo viene brevemente illustrato nei suoi aspetti fondamentali il funzionamento dell’interfaccia grafica da noi sviluppata in ambiente MATLAB e utilizzata nel seguito del capitolo per la simulazione dell’evoluzione temporale del modello di Rosenzweig-MacArthur. In particolare rivolgiamo la nostra attenzione all’interfaccia RM [7] che è stata creata specificatamente per il modello considerato e prevede la funzionalità aggiuntiva della valutazione dei punti fissi secondo i risultati analitici da noi riportati nel paragrafo 1.1.1, mentre segnaliamo l’esistenza dell’interfaccia avanzata DSSim [8] che costituisce una generalizzazione della precedente consentendo la simulazione di un qualsiasi sistema dinamico in tempo continuo; ovviamente la versatilità della seconda soluzione implica in generale un aumento, talvolta considerevole, del tempo necessario per il compimento della simulazione nonché l’assenza della funzionalità di calcolo dei punti fissi. Per maggiori dettagli, qui omessi, sul funzionamento delle due interfacce grafiche da noi realizzate, si faccia riferimento ai rispettivi manuali d’uso [7, 8] e alla documentazione di MATLAB sull’argomento [9, 10].

L’idea di realizzare l’interfaccia grafica RM è nata per verificare, in maniera facile e veloce, alcuni risultati ottenuti con CONTENT [11, 12], un tool molto potente per l’analisi dei sistemi dinamici in tempo continuo e tempo discreto (ODE, PDE, mappe iterate, DAE) che consente di effettuarne lo studio di biforcazione, ma che talvolta risulta poco user-friendly; è stata brevemente sperimentata anche l’interfaccia grafica di CONTENT per MATLAB (MATCONT [13]), sviluppata solamente per le equazioni differenziali ordinarie, che però non ha dimostrato gli indiscutibili pregi ed efficienza della versione classica indipendente.

L’interfaccia grafica RM consente di gestire in maniera immediata i parametri del modello dinamico, nonché quelli relativi alla simulazione del sistema che è stato definito mediante Simulink; infatti, la finestra principale dell’interfaccia grafica si occupa dell’impostazione dei parametri e della visualizzazione automatica dei risultati della simulazione, mentre un modello ausiliario di Simulink definisce il sistema dinamico attraverso una descrizione schematica a blocchi. Per quanto riguarda, invece, la valutazione dei punti fissi del sistema è stato previsto un pannello apposito che riceve i valori dei parametri del modello e provvede quindi alla visualizzazione delle coordinate numeriche dei punti fissi e dei relativi autovalori, calcolati mediante la matrice jacobiana associata al sistema; è possibile inoltre impostare automaticamente le condizioni iniziali del sistema sul pannello principale ai valori delle coordinate dei punti fissi semplicemente selezionando la riga di interesse nella finestra ausiliaria.


13

2.2 Simulazione dell’evoluzione temporale del modello di Rosenzweig-MacArthur: verifica delle tipologie di comportamento del sistema dinamico

In questo paragrafo riportiamo i risultati di numerose simulazioni condotte mediante l’interfaccia grafica RM circa l’evoluzione temporale del modello di Rosenzweig-MacArthur; in particolare abbiamo seguito il criterio, già illustrato nel paragrafo 1.2, di variare solamente i valori dei parametri r e K per verificare le diverse tipologie di comportamento del sistema dinamico individuate in [1] e brevemente considerate in [6]. 2.2.1 Estinzione del superpredatore Come abbiamo già visto nel paragrafo 1.2, esiste una vasta zona nello spazio ( K , r ) per la quale si assiste all’estinzione del superpredatore. Occorre subito evidenziare che, soprattutto in questo caso, la tipologia di comportamento del sistema non è considerata tanto nel suo aspetto qualitativo complessivo, quanto nella sua significatività a livello di dinamica di popolazione; si può infatti notare che sono state considerate in blocco zone dello spazio dei parametri accomunate dall’assenza di attrattori o cicli stabili che prevedano la sopravvivenza della popolazione 3x , mentre la tipologia delle traiettorie varia sensibilmente in funzione dei valori dei parametri r e K . Nella prima figura riportiamo il risultato della simulazione ottenuta utilizzando l’impostazione dei parametri di default dell’interfaccia RM [9]; come si può notare, il superpredatore si estingue, mentre le concentrazioni di preda e predatore oscillano in modo regolare garantendone la sopravvivenza.


14

Riportiamo ora i risultati della valutazione dei punti fissi: in particolare esistono due punti sella (1 e 2), peraltro sempre presenti (vedi paragrafo 1.1.1), e un fuoco instabile (3) che tende ad attrarre ed espandere le traiettorie sul piano

03 =x ; si noti che le tre soluzioni per i punti fissi sono valide anche per un eventuale sistema ridotto che modelli una catena alimentare costituita solamente da una preda ed un predatore.

Il risultato riportato nelle due figure seguenti è stata ottenuto impostando le condizioni iniziali ai valori delle coordinate del punto fisso 3; è evidente che l’impossibilità di esprimere le coordinate del punto fisso con precisione infinita, determina la divergenza della traiettoria dal fuoco instabile verso un regime ciclico di coesistenza di preda e predatore.


15

Come abbiamo già avuto modo di evidenziare, l’estinzione del superpredatore può avvenire in diverse condizioni di coesistenza tra preda e predatore; nel caso seguente, infatti, si può notare, oltre alla presenza dei due punti sella 1 e 2, l’esistenza del fuoco stabile 3.


16

L’esempio seguente è stato realizzato per verificare i risultati della ricerca analitica dei punti fissi riportati nel paragrafo 1.1.1; attraverso calcoli più o meno agevoli abbiamo potuto individuare, per 1=K e 53.047.0 << r (i valori limite di r sono approssimati), una ristretta zona dei parametri con cinque soluzioni per i punti fissi.


17

Riportiamo per completezza le soluzioni dei punti fissi in corrispondenza delle impostazioni dei parametri utilizzate per la precedente simulazione; la soluzione 4- evidenziata è stata utilizzata come condizione iniziale.

2.2.2 Coesistenza stazionaria Nella figura seguente è riportato il risultato di una simulazione effettuata con valori dei parametri e delle condizioni iniziali standard per la verifica del comportamento stazionario del sistema.


18

Accanto ai consueti punti sella (1 e 2), verifichiamo la presenza di un fuoco instabile (3) e di un punto fisso non iperbolico (4+) con due autovalori complessi a parte reale negativa ed uno nullo; come è noto dalla teoria [14], la presenza di un’autovalore nullo per il sistema linearizzato non consente di valutare a priori la natura del punto fisso, ma i risultati della simulazione precedente e di quella seguente manifestano un comportamento globalmente attrattivo (le condizioni iniziali dell’esempio seguente sono state impostate come piccola perturbazione su 3x delle coordinate del fuoco instabile 3).


19

Nella figura seguente è riportato il risultato di una simulazione condotta annullando la condizione iniziale per il superpredatore; con non troppa sorpresa abbiamo ottenuto un comportamento all’apparenza stabile che è valido però nel caso di due sole popolazioni in competizione, dal momento che nel sistema ridotto di ordine 2 scompare l’autovalore reale positivo: infatti, nel nostro caso, basterebbe una condizione iniziale per 3x leggermente diversa da zero per ricondurci al comportamento standard dei due esempi precedentemente illustrati.

Per quanto riguarda la coesistenza stazionaria riportiamo un ultimo esempio, in cui una concentrazione iniziale di 3x nettamente superiore rispetto a quella delle altre due specie determina, in prima istanza, una sensibile diminuzione delle concentrazioni 1x e 2x ; non appena la densità del superpredatore diminuisce, la concentrazione di preda e predatore può aumentare fino al raggiungimento dell’equilibrio 4+ già individuato.


20

2.2.3 Coesistenza ciclica in bassa frequenza I risultati delle simulazioni di questo paragrafo sono stati ottenuti fissando i parametri r e K e variando opportunamente le condizioni iniziali; nella figura seguente è possibile notare la tipica forma di un teacup attractor, già descritta nel paragrafo 1.2, dovuta alla contemporanea presenza di oscillazioni a bassa frequenza per la concentrazione del superpredatore e ad alta frequenza per le densità di preda e predatore [1].


21

Le due figure seguenti sono state ottenute impostando le condizioni iniziali delle popolazioni rispettivamente alle coordinate dei punti fissi 3 e 4+; nel primo caso la struttura dell’attrattore viene completamente alterata dalla condizione iniziale nulla per 3x , mentre nel secondo caso si può individuare una forma abbastanza simile al tipico teacup attractor.


22

2.2.4 Coesistenza caotica Anche nel caso della coesistenza caotica, tutte le simulazioni di cui riportiamo i risultati sono state effettuate mantenendo fissati i valori di r e K e variando di volta in volta le condizioni iniziali.


23

Nelle due figure precedenti e nel dettaglio qui sotto è possibile individuare la forma classica di un cut-teacup attractor che evolve con comportamento caotico in una regione limitata; in particolare l’evoluzione caotica è più facilmente valutabile analizzando i risultati delle evoluzione temporale delle concentrazioni delle tre popolazioni, in particolare quella del superpredatore.

Abbiamo quindi attivato il calcolo dei punti fissi per la configurazione dei parametri già utilizzata nelle due precedenti simulazioni; la simulazione del sistema con condizioni iniziali uguali alle coordinate del punto fisso 4+ ha dato risultati analoghi ai due casi precedenti (vedi figura pagina successiva).


24

2.2.5 Coesistenza ciclica in alta frequenza Per illustrare il comportamento ciclico in alta frequenza riportiamo due simulazioni; nel primo caso la condizione iniziale per la concentrazione del superpredatore è superiore ai valori della coesistenza ciclica e pertanto l’assestamento del comportamento è visibile solamente alla fine di un periodo di transitorio.


25

La seconda simulazione è stata effettuata variando sia il valore di r che la condizione iniziale per 3x e, come nel caso precedente, la coesistenza ciclica si manifesta solamente alla fine di un periodo di transitorio; in entrambi i casi l’attrattore ha perso la forma del cut-teacup attractor.


26

3.0 ANALISI DI BIFORCAZIONE Nel capitolo precedente abbiamo verificato, attraverso la simulazione temporale del sistema, la presenza delle diverse tipologie di comportamento del modello di Rosenzweig-MacArthur individuate in [1, 2]; per fare ciò è stato necessario fissare di volta in volta il valori dei parametri r e K , mentre gli altri otto parametri sono sempre stati considerati fissati ai valori indicati nel paragrafo 1.2 . Vogliamo ora dedicare la nostra attenzione all’analisi delle biforcazioni locali del sistema cercando di ottenere mediante CONTENT gli stessi risultati già ottenuti con AUTO [5] in [1, 2]; anche in questo caso le nostre simulazioni sono state effettuate variando i valori di r e K e mantenendo fissati i restanti parametri.


27

3.1 Generalità sull’utilizzo di Content CONTENT [11, 12] è un tool software (ne esistono varie versioni per diverse piattaforme) orientato all’analisi dei sistemi dinamici e delle loro proprietà e dotato di un insieme di apposite funzionalità che consentono di effettuare lo studio di biforcazione dei punti fissi al variare dei parametri. Il programma, fornito di un interfaccia grafica per la scelta delle impostazioni e dei parametri di simulazione e la successiva rappresentazione grafica dei risultati, consente di analizzare sistemi dinamici descritti da modelli matematici a tempo continuo (equazioni differenziali ordinarie, equazioni alle derivate parziali) e a tempo discreto (mappe iterate) mediante la loro espressione analitica; dopo aver provveduto alla definizione del sistema è possibile effettuare diverse tipologie di simulazione in funzione delle necessità: evoluzione temporale del sistema, ricerca dei punti fissi e valutazione della loro stabilità, diagrammi di biforcazione. Esistono poi strumenti ausiliari che consentono una più approfondita valutazione dei risultati delle simulazioni, come ad esempio il pannello numerico, che visualizza i valori istantanei di variabili di stato, parametri ed autovalori, o le varie funzionalità grafiche per la gestione e la rappresentazione dei risultati. CONTENT è indiscutibilmente uno strumento molto potente, ma non si presta agevolmente alla realizzazione di numerose simulazioni consecutive ed all’esportazione dei risultati grafici delle stesse; per questi motivi le simulazioni temporali ampiamente descritte e documentate nel precedente capitolo sono state opportunamente verificate mediante CONTENT, ma i risultati riportati sono quelli ottenuti utilizzando l’interfaccia grafica da noi sviluppata in ambiente Matlab. Inoltre, come vedremo, abbiamo ricavato grandi vantaggi dall’utilizzo della nostra funzionalità di calcolo degli autovalori del sistema linearizzato ai fini dell’individuazione di regioni interessanti nello spazio dei parametri, cioè quelle zone in cui il numero dei punti fissi è più consistente e i diagrammi di biforcazione risultano conseguentemente più complessi e ricchi di informazioni circa il comportamento del sistema dinamico nel suo complesso. Nel paragrafo seguente descriveremo quindi il nostro approccio allo studio di biforcazione del modello di Rosenzweig-MacArthur illustrandone i comportamenti generali e confrontandolo con quanto abbiamo già riportato nel paragrafo 1.2 circa la struttura di questo sistema dinamico che è già stato ampiamente analizzato e descritto in letteratura.


28

3.2 Risultati dell’analisi In questo paragrafo presentiamo i risultati dell’analisi di biforcazione per il modello di Rosenzweig-MacArthur condotta mediante CONTENT; come già anticipato, abbiamo scelto di analizzare una regione dello spazio (r , K ) interessata dalla presenza di soluzioni valide per tutti i punti fissi ricavati per via analitica (paragrafo 1.1.1) e ci siamo concentrati sui risultati ottenuti dalle simulazioni già illustrate nel precedente capitolo circa l’utilizzo dell’interfaccia MATLAB da noi realizzata. Inizialmente abbiamo deciso di considerare i risultati della valutazione dei punti fissi ottenuti per 5.0=r e 0.1=K (ultima simulazione del paragrafo 2.2.1) che indicano la presenza di un punto fisso di coordinate ( )5.7,125.0,5.0 : abbiamo così potuto impostare tutti i dati per la ricerca dei punti di equilibrio a partire dal punto fisso considerato per K fissato e r variabile.

La figura precedente evidenzia sette punti “particolari” etichettati in funzione delle loro proprietà; attraverso un’analisi più approfondita, condotta mediante l’utilizzo dei risultati numerici e dei diagrammi di biforcazione delle concentrazioni delle tre popolazioni in funzione dei parametri r e K , ci occuperemo di verificarne la consistenza e la coerenza con i risultati già a nostra disposizione.


29

Dai sei grafici precedenti si può subito notare che tre dei sette punti evidenziati nel diagramma di biforcazione sul piano (r , K ) non sono interessanti ai fini della nostra analisi, dal momento che almeno una delle loro coordinate risulta negativa. Inoltre, a conferma di ciò, abbiamo verificato che, per valori dei parametri compresi nell’intervallo di nostro interesse ( 20 ≤≤ r , 20 ≤< K ), almeno uno dei rami dei diagrammi di biforcazione ottenuti da ciascuno dei tre suddetti punti permane fuori dell’ottante positivo delle concentrazioni delle tre popolazioni: in particolare il punto di biforcazione P6 coincide con il quinto punto fisso ricavato per via analitica nel paragrafo 1.1.1 e immediatamente scartato; i tre punti (contraddistinti con P4, P5 e P6 nella figura della pagina precedente) non saranno pertanto più considerati nel seguito della nostra analisi.


30

Per una migliore comprensione dei grafici precedenti abbiamo riunito nella seguente tabella tutti i dati relativi ai punti evidenziati da CONTENT; per ciascun punto sono riportati il numero distintivo da noi assegnato nella figura di pag. 28, la tipologia della biforcazione, i valori del parametro r (in questa prima simulazione K è stato mantenuto fisso a 1) e le coordinate dei punti fissi in corrispondenza della stessa biforcazione: si noti che, per tutti i punti del diagramma, la concentrazione del predatore è costante e pari a 125.0 , cioè il valore di Y del paragrafo 1.1.1.

N. Tipo r X1 X2 X3 Note P.F.

MatlabP1 LP 0.46875 0.3332954 0.125 5.416074 Comparsa punti fissi 4+ e 4- 4+ 4-P2 H 0.5200243 0.5426707 0.125 7.905922 Inizio regione stazionaria 4+P3 BP 0.5308824 0.1052631 0.125 0 Scomparsa punto fisso 4- 3P4 H 0.5765284 0.004508673 0.125 -2.51782 Punto fisso negativoP5 H 0.6236658 0.001071383 0.125 -4.933253 Punto fisso negativoP6 BP 0.6250001 0 0.125 -5.000004 Punto fisso negativoP7 H 0.7612113 0.7465619 0.125 9.402668 Inizio regione ciclica 4+

Se si osservano i dati contenuti nella tabella precedente e i sei grafici di pagina 29, si possono chiaramente comprendere le condizioni che determinano la presenza dei punti fissi 4+ e 4-, già ampiamente analizzati nel paragrafo 1.1.1 e qui valutati per 0.1=K :

§ per ( ) 333

332

2

24

deadb

bK

Kar

−+< i punti fissi 4+ e 4- non esistono;

§ per ( )

46875.04

333

332

2

2 ≅−+

=dea

db

bK

Kar i punti fissi 4+ e 4- coincidono e il punto

fisso ha le coordinate del punto limite P1 indicate in tabella;

§ il punto fisso 4+ esiste sempre per ( ) 333

332

2

24

deadb

bK

Kar

−+> tendendo alle

coordinate asintotiche

−

+=

−== 2

2

22

3

3

333

33

1,, d

Xbea

Yde

Zdea

dbYKX che

nel caso specifico valgono ( )625.10,125.0,1 ;

§ il punto fisso 4- esiste per ( ) 333

332

2

24

deadb

bK

Kar

−+> fino a quando è soddisfatta

la condizione Hdea

dbX =

−≥

222

22 (l’uguaglianza è soddisfatta in

corrispondenza del punto P3 dove i punti 3 e 4- coincidono). Ci dedichiamo ora alla continuazione dei diagrammi di biforcazione a partire dai punti di Hopf P2 e P7; in particolare ci concentriamo sul punto P2 che viene selezionato come punto iniziale di simulazione per r e K variabili, anche se lo stesso risultato sarebbe raggiunto partendo da P7.


31

Se confrontiamo il grafico appena ottenuto con il diagramma di riferimento per le diverse tipologie di comportamento del modello di Rosenzweig-MacArthur che abbiamo riportato a pagina 9, scopriamo che i due punti di Hopf P2 e P7 determinano il passaggio rispettivamente tra la condizione di morte del superpredatore e la coesistenza stazionaria delle tre specie e tra la coesistenza stazionaria e il regime ciclico in bassa frequenza [1]; per questo motivo la curva che definisce lo spostamento di tali punti nel piano (r , K ) rappresenta le linee di confine tra le suddette regioni. Nel grafico precedente CONTENT segnala la presenza di tre punti di “biforcazione generalizzata di Hopf”, cioè punti di Hopf per i quali si ha una coppia di autovalori complessi a parte reale nulla (si veda a proposito l’help del programma); su tali punti non concentreremo la nostra attenzione, dal momento che il loro studio comporterebbe la necessità di considerare la variazione di un terzo parametro. Nei sei grafici seguenti sono mostrati, oltre alle curve originarie della prima simulazione, gli andamenti delle coordinate dei punti di Hopf al variare dei parametri r e K ; si noti in particolare la concentrazione costante del predatore in corrispondenza dei suddetti punti.


32

Concentriamo ora la nostra attenzione allo studio della variazione delle coordinate del punto limite P1 in funzione dei parametri r e K ; come abbiamo già detto, il punto P1 determina la simultanea nascita dei punti fissi 4+ e 4-, mentre mediante altre considerazioni circa il segno dei risultati si determina la loro significatività fisica per il sistema dinamico che stiamo considerando. CONTENT segnala inoltre la presenza di tre punti particolari: due “cuspidi”, cioè punti con un autovalore nullo e tutti gli altri con parte reale diversa da zero, e una biforcazione di Bogdanov-Takens, cioè punto avente due autovalori nulli;


33

anche in questo caso, per maggiori dettagli sulla classificazione dei punti, rimandiamo all’help del programma.

Nella figura qui sopra e nelle sei successive sono riportate, sugli stessi grafici che contengono i risultati delle precedenti analisi, le curve che descrivono l’andamento delle coordinate del punto limite P1; a questo proposito notiamo che la concentrazione del predatore si mantiene costante per qualsiasi valore di r e K e soprattutto che lo spostamento del punto limite può determinare un radicale cambiamento della distribuzione dei punti fissi, come la completa scomparsa di una soluzione valida per il punto 4- o la variazione dell’intervallo dei parametri per cui esiste la soluzione 4+.


34

Abbiamo già visto (pag. 30) che i punti fissi 3 e 4- sono coincidenti per particolari configurazioni dei parametri del modello; ciò è confermato dai risultati di simulazione che riportiamo di seguito: in particolare si noti che la curva che descrive le coordinate del punto fisso in funzione di r segue l’andamento prevedibile dalle espressioni ricavate nel paragrafo 1.1.1.


35


36

4.0 CONCLUSIONI L’analisi del modello di Rosenzweig-MacArthur da noi condotta si è mossa secondo tre linee guida principali: ricerca analitica dei punti fissi, verifica delle tipologie di comportamento mediante MATLAB, analisi di biforcazione con CONTENT; in particolare abbiamo cercato di verificare la coerenza delle indicazioni ottenute dai diversi approcci, ottenendo risultati soddisfacenti per la definizione della struttura complessiva del sistema dinamico. Abbiamo potuto scoprire “sul campo” le proprietà tipiche di sistemi dinamici complessi che manifestano molteplici tipologie di comportamento, la compresenza di attrattori, la comparsa del caos ed un’estrema sensibilità alle condizioni iniziali, e abbiamo verificato che la coesistenza delle tre popolazioni tende a variare la sua complessità in funzione della capacità di autosostentamento della preda andando a massimizzare il rendimento delle risorse in corrispondenza del passaggio dal comportamento caotico a quello ciclico in alta frequenza [1, 2]. In conclusione vogliamo ricordare che nella nostra analisi ci siamo limitati a studiare l’impatto sul sistema dinamico di solamente due dei dieci parametri [1, 2, 3] del modello di Rosenzweig-MacArthur che potrebbe pertanto descrivere situazioni interessanti anche per differenti configurazioni, non solo dal punto di vista puramente teorico, ma anche a livello biologico; segnaliamo inoltre la possibilità di effettuare lievi modifiche al modello [2, 15] per considerare effetti minimizzati od ignorati nella descrizione del sistema che abbiamo considerato.


37

Riferimenti bibliografici

[1] A. Gragnani, O. De Feo, S. Rinaldi, FOOD CHAINS IN THE CHEMOSTAT: RELATIONSHIPS BETWEEN MEAN YIELD AND COMPLEX DYNAMICS, Bullettin of Mathematical Biology, vol. 60, pp. 703-719, 1998 http://ideas.repec.org/p/wop/iasawp/ir97042.html

[2] O. De Feo, S. Rinaldi, YIELD AND DYNAMICS OF TRITROPHIC FOOD CHAINS, The American Naturalist, vol. 150 n° 3, pp. 328-345, 1997 http://www.journals.uchicago.edu/AN/journal/contents/v150n3.html

[3] Yu. A. Kuznetsov, S. Rinaldi, REMARKS ON FOOD CHAIN DYNAMICS, Mathematical Biosciences, vol. 134, pp. 1-33, 1996 http://www.sciencedirect.com/science/journal/00255564

[4] M.L. Rosenzweig, PARADOX OF ENRICHMENT: DESTABILIZATION OF EXPLOITATION

ECOSYSTEMS IN ECOLOGICAL TIME,. Science, vol. 171, pp. 385-387, 1971 [5] E. Doedel, J. Kernévez, AUTO: SOFTWARE FOR CONTINUATION PROBLEMS IN

ORDINARY DIFFERENTIAL EQUATIONS WITH APPLICATIONS, California Institute of Technology Technical Report, Applied Mathematics, CALTEC, Pasadena, 1986 http://www.enm.bris .ac.uk/staff/hinke/dss/continuation/auto.html

[6] A. Basteris, M.Izzo, Tesi di laurea triennale: DINAMICHE COMPLESSE IN SISTEMI

PREDATORE, Università degli Studi di Genova, Novembre 2003 [7] E. Costa, P. Spiga, MANUALE D’USO PER “INTERFACCIA GRAFICA IN AMBIENTE

MATLAB PER LA SIMULAZIONE DELL'EVOLUZIONE TEMPORALE E LA VALUTAZIONE DEI

PUNTI FISSI DI UN SISTEMA DINAMICO NON LINEARE RAPPRESENTANTE IL MODELLO

PREDA - PREDATORE - SUPERPREDATORE DI ROSENZWEIG -MACARTHUR”, Genova, 2004

[8] E. Costa, P. Spiga, MANUALE D’USO PER “INTERFACCIA GRAFICA IN AMBIENTE

MATLAB PER LA SIMULAZIONE DELL'EVOLUZIONE TEMPORALE DI UN SISTEMA

DINAMICO DEFINITO MEDIANTE S-FUNCTION IN UN MODELLO SIMULINK”, Genova, 2004

[9] MATLAB website: http://www.mathworks.com/

[10] The MathWorks, CREATING GRAPHICAL USER INTERFACES – VERSION 6, 2000-2002 http://www.mathworks.com/access/helpdesk/help/pdf_doc/matlab/buildgui.pdf

[11] CONTENT website http://www.math.uu.nl/people/kuznet/CONTENT/

[12] Yu. A. Kuznetsov, CONTENT – INTEGRATED ENVIRONMENT FOR ANALYSIS OF

DYNAMICAL SYSTEMS. TUTORIAL, 1998 http://www.enm.bris.ac.uk/staff/hinke/dss/continuation/content.html

[13] MATCONT website http://allserv.rug.ac.be/~ajdhooge/research.html

[14] S.H. Strogatz, NONLINEAR DYNAMICS AND CHAOS WITH APPLICATIONS TO PHYSICS ,

BIOLOGY, CHEMISTRY AND ENGINEERING, Addison-Wesley Publishing Company, 1994 [15] O. De Feo, S. Rinaldi, SINGULAR HOMOCLINIC BIFURCATIONS IN TRI-TROPHIC FOOD

CHAINS, Mathematical Biosciences, vol. 148, pp. 7-20, 1997 http://www.sciencedirect.com/science/journal/00255564

Corso di Sistemi caotici per la protezione dell’informazione...

Documents

Transcript of Corso di Sistemi caotici per la protezione dell’informazione...