Modellistica e Simulazione: il sistema di Rosslera.a. 2007-2008 Chiara Mocenni

Il sistema di Rossler è il più semplice sistema di terzo ordine capace

di manifestare comportamenti di tipo caotico.

( )

x y zy x ayz b z x c

= − −= += + −

Questo sistema presenta un solo termine non lineare, il prodotto tra z

e x nella terza equazione, il quale è parte fondamentale del processo di

stretching and folding tipico dei sistemi caotici.

Consideriamo le prime due equazioni e supponiamo che z sia

abbastanza piccolo da essere trascurato. Pertanto il sottosistema

x yy x a= −= + y

(1)

può essere trasformato nell’oscillatore lineare del secondo ordine

( )22

1,2

0

4

2

x ax x

a aλ

− + =

± −=

Con il parametro a positivo questo oscillatore ha uno smorzamento

negativo, in più se a è compreso tra 0 e 2 l’origine è una spirale instabile.

Vedendo il tutto in tre dimensioni questo significa che traiettorie vicine al

piano (x,y) saranno respinte in maniera spiraliforme dall’origine. Questo

genera il dispiegamento delle traiettorie adiacenti, il quale è il primo

ingrediente del caos

Tutto ciò è generato solo da termini lineari, ma questo non basta per

ottenere un comportamento caotico perché, se l’intero sistema fosse

lineare, l’effetto di stretching continuerebbe fino a far divergere le

traiettorie all’infinito. Per confinare l’azione di allontanamento all’interno

di una regione di spazio limitata, il termine non lineare è necessario.

Il parametro c nella terza equazione agisce da soglia per attivare

l’azione non lineare di ripiegamento o folding. Considerando solamente la

terza equazione si nota che, quando il valore di x è minore della costante c,

il coefficiente di z è negativo, il sottosistema z è stabile e tende a riportare

z a valori prossimi a

( )b x c− −

-b/(x-c)

dz

z

b

D’altro canto, se x supera c allora z appare nella terza equazione

moltiplicato per un fattore positivo e il precedente sistema stabile diverge.

Scegliendo b>0 la divergenza avviene nel senso positivo degli z come è

possibile vedere nella figura seguente.

Cercheremo a questo punto di descrivere nella maniera più completa

possibile i comportamenti del sistema in funzione del parametro più

significativo e la scelta è ricaduta sul parametro a poiché la sua variazione

determina la dinamica del sottosistema (1), cioè per certi dei suoi valori

l’oscillatore di secondo ordine presenta un punto fisso stabile oppure

instabile. Cercheremo a questo punto di capire come tutto ciò si rifletta sul

sistema considerato nella sua interezza e lo faremo fissando b e c

rispettivamente ai valori 2 e 4 e calcolando i punti fissi del sistema in

funzione del parametro a. I valori di b e di c sono stati scelti poiché

rendono semplici molti conti, comunque, dove possibile, sarà specificato

se lo studio fatto vale in generale o solo per i particolari valori di b e c.

( )( )

( )

21,2

21,2

21,2

1 42

1 42

1 42

x c c ab

y c c aba

z c c aba

= ± −

= − ± −

= ± −

Per b e c positivi, il sistema presenta due punti fissi fintanto che a

è diverso da zero e minore di c2/(4b).

E’ utile quindi concentrasi inizialmente su valori di a negativi, vedere poi

cosa succede se a=0 per poi spostarsi sull’intervallo 0<a< c2/(4b).

Per fare ciò è necessario calcolare la matrice A del sistema

linearizzato che risulta uguale a:

0 1 11 0

0A a

z x

− −

c

⎡ ⎤⎢ ⎥= ⎢ ⎥⎢ ⎥−⎣ ⎦

Fissato un particolare valore di a negativo, ad esempio -6, e

utilizzando un piccolo script su MATLAB, sono stati calcolate le

coordinate dei punti, gli autovalori della matrice del sistema linearizzato A

e gli autovettori associati. Dai risultati ottenuti si nota che uno dei due

punti fissi è un fuoco stabile mentre l’altro una sella come si può vedere

dai dati e dalla figura seguenti. Da ora in poi chiameremo la sella S1ed il

punto fisso P2, tanto che abbia un comportamento stabile che instabile.

a=-6

b=2

c=4

z1=(0.5/a)*(c+sqrt(c^2-4*a*b));

z2=(0.5/a)*(c-sqrt(c^2-4*a*b));

y1=-z1;

y2=-z2;

x1=a*z1;

x2=a*z2;

punto_fisso_1=[x1 y1 z1]

A1=[0 -1 -1;1 a 0;z1 0 x1-c]

[AUTOVETTORI_1,DELTA_1]=eig(A1)

punto_fisso_2=[x2 y2 z2]

A2=[0 -1 -1;1 a 0;z2 0 x2-c]

[AUTOVETTORI_2,DELTA_2]=eig(A2)

punto_fisso_1 = 6 1 -1

A1 = 0 -1 -1

1 -6 0

-1 0 2

AUTOVETTORI_1 = -0.3689 -0.9186 0.1729

-0.0439 -0.1693 0.9847

0.9284 -0.3570 0.0221

AUTOVALORI_1 = 2.3974 0 0

0 -0.5729 0

0 0 -5.8245

punto_fisso_2 = -2.0000 -0.3333 0.3333

A2 = 0 -1.0000 -1.0000

1.0000 -6.0000 0

0.3333 0 -6.0000

AUTOVETTORI_2 = 0.9837 -0.2142 0.0000

0.1705 -0.9267 -0.7071

0.0568 -0.3089 0.7071

AUTOVALORI_2 = -0.2311 0 0

0 -5.7689 0

0 0 -6.0000

Sella S1

x=6 y=1 z=-1

Fuoco Stabile P1

x=-2 y=-0.3333 z=0.33333

Figura 1 : Fuoco stabile e sella per a=-6 b=2 c=4

Successivamente i valori dei due punti sono stati inseriti in

MATCONT con il quale è stato fatto lo studio di biforcazione per valori di

a minori di zero. Per a che tende a meno infinito non sono state riscontrate

variazioni nel comportamento qualitativo dei due punti, l’unico tipo di

differenza si nota nella dislocazione spaziale di S1 e P2 poiché in entrambi

si ha convergenza a zero lungo gli assi y, z e una divergenza verso infinito

lungo l’asse x. Il punto H riportato da MATCONT sulla traiettoria di

biforcazione del punto di sella rappresenta una ‘neutral saddle’, cioè S1

per a=-2.69929 presenta due autovalori reali uguali in modulo ma di segno

opposto.

Sella neutra 2 autovalori uguali ma di segno opposto

Sella S1

Fuoco Stabile P2

Entrambi i punti divergono verso infinito lungo x

Figura 2 : Diagramma di biforcazione per a minore di zero (x contro a)

Figura 3 : Diag

Per a=0 si ha un

poiché i punti fissi non

mentre il punto stabile

all’infinito. Per verificar

due punti.

1

1

2

2

2

l

l

l

a

a

z

z

z

y

x

=

=

=

=

= −

Per a -> 0 P2 converge ad un valore finito mentre S1 diverge all’infinito lungo y e z e converge a

Per a -> -∞ i due punti convergono a zero lungo y e z e divergono lungo x

ramma di biforcazione per a minore di zero (y contro a)

comportamento del sistema abbastanza bizzarro

mutano il loro comportamento qualitativo ma,

P2 presenta coordinate finite, la sella S1 diverge

e ciò basta fare il limite per a che tende a zero dei

( )( )( )

2

0

2

0

2

0

2

2

1im 421im 4

21im 4

2

0

a

c c aba

c c aba

bc c aba c

bzc

az

−

+

→

→

→

+ − =

+ − =

− − =

− = −

=

+∞

−∞

Di seguito a tutto ciò, poiché il sistema presenta un solo punto fisso

stabile, si è cercato di vedere se la stabilità fosse di tipo globale o solo

locale. Purtroppo, da risultati sperimentali, si vede bene che esistono tre

bacini di attrazione, i quali sono uno associato al punto stabile mentre gli

altri ad una struttura a questo punto più complessa di un punto fisso che fa

divergere le traiettorie all’infinito.

3° bacino di attrazione

2° bacino di attrazione

1° bacino di attrazione

Per a>0 si ha la dinamica più interessante dato che il sistema si

avvicina alla zona di comportamento caotico. Comunque inizialmente

prenderemo in considerazione il punto S1, il quale ricompare dall’infinito

e continua a comportarsi come una sella fin tanto che a risulta minore di

c2/(4b). Al raggiungimento di questo valore si ha una biforcazione di tipo

nodo-sella, cioè S1 e P2 collidono in un unico punto per poi scomparire

definitivamente per a> c2/(4b). Questo si vede bene tanto dalle simulazioni

con MATCONT che dallo studio dal punto di vista analitico dei punti fissi,

infatti quando a=c2/(4b) si ha che:

22

1,2

1,2

1,2

4 04

12

1 22

1 22

ca c abb

x c

by ca c

bz ca c

= → − =

=

= − = −

= =

Punto P2

Sella S1

Biforcazione Nodo-Sella

Nella figura precedente si può notare come MATCONT segnali il

punto di collisione con la sigla LP (Limit Point) per indicare che, superato

quel valore particolare di a, il punto fisso in esame scompare. Sovrapposto

al all’indicatore LP è presente l’ulteriore marker H, il quale indica che il

punto fisso, risultante dalla fusione di S1 e P1 in a= c2/(4b), è una “sella

neutra” (‘Neutral Saddle’), cioè ha due autovalori uguali in modulo ma di

segno opposto.

A questo punto è necessario concentrarsi sulle biforcazioni cui va in

contro il punto P2 per a compreso nell’intervallo [0,c2/4b] nel caso

particolare di b=2 e c=4. Dallo studio fatto con MATCONT si è visto che

il punto presenta una biforcazione di Hopf supercritica quando

a=0.12496748. Conseguentemente P2 diviene una spirale instabile (due

autovalori complessi e coniugati attraversano lo zero divenendo a parte

reale positiva mentre il terzo autovalori è reale negativo) circondata un

ciclo limite stabile.

Fuoco Stabile

Biforcazione di Hopf

Spirale Stabile

Figura 4 : Biforcazione do Hop

Dopo la biforcazione di Hopf, il p

precedenza, è divenuto una spirale insta

biforcazione nodo-sella.

A questo punto le coordinate del punt

MATCONT il quale ci ha permesso di trov

ciclo, il suo periodo e di provarne effettivam

Sella

Biforcazione Nodo-Sella

Spirale Instabile

f (x contro a)

unto P2, come gia detto in

bile e così rimane sino alla

o di Hopf sono state inserite in

are la collocazione spaziale del

ente la stabilità.

Spirale Instabile T=6.19324 X=0.1026 Y=-0.5132 Z=0.5132 a=0.2

Ciclo Limite

Figura 5 : Ciclo limite per a=0.2

Come si può vedere dalla figura precedente la spirale instabile P2 è

circondate da un nuovo attrattore, il ciclo limite C3, il quale è stabile e in

questo caso, cioè con a=0.2, ha periodo T=6.19324. Per a proprio uguale al

valore di biforcazione il periodo T del ciclo è 6.233078.

Ci concentriamo adesso sulle possibili evoluzioni del ciclo limite C3

e per fare ciò utilizziamo sempre MATCONT, posizionandoci sulle

coordinate del punto di Hopf e continuando a far variare il parametro a.

Poiché ci avviciniamo sempre di più alla zona di comportamento caotico ci

aspettiamo che il ciclo limite presenti dei branching o dei raddoppiamenti

di periodo che lo spingano verso il caos.

Evoluzione del ciclo limite al variare di a

Raddoppiamento di periodo

Figura 6 :Raddoppiamento di periodo per a=0.3348

Come si vede dalla figura e come ci aspettavamo, il ciclo

all’aumentare di a presenta un raddoppiamento di periodo per a=0.3348.

Si comprende bene come questo primo raddoppiamento sia legato in

qualche modo all’azione della terza equazione del sistema, questo perché

all’aumentare del parametro a, il ciclo limite C3, che è più o meno centrato

in zero, aumenta di raggio ed in particolare vicino al valore del

raddoppiamento di periodo il ciclo limite raggiunge il valore 4 lungo la

coordinata x. Questo valore, per c=4, destabilizza la terza equazione del

sistema mettendo in moto il processo di folding gia citato nelle prime

pagine e causando la divisione dell’orbita.

Ciclo Limite per a=0.2

Ciclo Limite per a=0.32

Figura 7 : All'aumentare del valore del parametro il ciclo limite si avvicina al valore di inversione di stabilità della terza equazione del sistema

Questo si comprende meglio pensando che l’orbita del ciclo limite

attraversa due zone ben distinte, separate dalla retta x=4 e di stabilità

opposta, cioè una stabile e una instabile lungo l’asse z. La penetrazione

dell’orbita nella zona instabile porta questa a divergere lungo le zeta

positive e continuerebbe all’infinito se non fosse presente una retroazione

negativa nelle prima equazione capace di riportala nella zona “di

equilibrio”.

Divergenza e successivo ripiegamento della traiettoria lungo z

Figura 8 : Ciclo di periodo raddoppiato per a=0.35

Comunque, ritornando all’analisi con MATCONT, si comprende

meglio ciò che succede se facciamo continuare il punto di raddoppiamento

di periodo. Il programma individua una orbita particolare che definisce

“branching point cycle”, superata la quale il ciclo limite C3 diviene

instabile. Come effetto collaterale si ha la formazione di un nuovo ciclo

limite stabile e di periodo doppio rispetto all’originale; in pratica il punto

di branching si comporta come una biforcazione di Hopf per i cicli.

Branchig Point Cycle

Evoluzione instabile del ciclo limite C3 dopo il Branching

Figura 9 : Andamento in funzione di a del ciclo limite C3 reso instabile dal Branching Point Cycle

Figura 10 : Andamento grafico x contro a

Branchig Point Cycle

Ciclo limite stabile di periodo doppio

Figura 11 : Evoluzione del ciclo limite di periodo doppio dopo il BPC

Ciclo limite instabile

Ciclo stabile di periodo doppio

Figura 12 : Coesistenza del ciclo limite instabile C3 e del suo raddoppiamento di periodo

Le figure precedenti mostrano l’andamento nel tempo del sistema nel

caso di stato iniziale prossimo al ciclo limite C3 e per un valore del

parametro a pari a 0.35. Si vede bene come la traiettoria inizialmente

rimane sul ciclo limite instabile per poi allontanarsi ed essere attratta verso

quello stabile e di periodo doppio. Questo comportamento mette in luce

nuovamente l’effetto di stretching e folding presente nel sistema, poiché

l’evoluzione stabile del ciclo C3 non è altro che un suo stiramento e

ripiegamento nello spazio.

Purtroppo dopo il primo raddoppiamento di periodo MATCONT non

ha più la sensibilità necessaria per detectare raddoppiamenti di periodo

successivi, pertanto da questo punto in poi l’analisi sarà portata avanti

basandosi solamente sui dati ottenuti dall’integrazione numerica.

Comunque, estrapolando il comportamento del sistema dai dati ottenuti

sino ad ora, ci aspettiamo una dinamica a cascata attraverso

raddoppiamenti di periodo successivi e sempre più frequenti al variare di a

e la conseguente generazione di cicli limite instabili. Pertanto esisterà un

valore di accumulazione per a, tale che, una volta sorpassato, il sistema

subentra nella zona di comportamento caotico.

Per individuare questo valore e circoscrivere la zona caotica, le

simulazioni sono state portate avanti in maniera ricorsiva, cioè attraverso

aggiustamenti successivi del parametro e conseguenti simulazioni. Il

grafico sottostante mostra i comportamenti del sistema relativi a vari

intervalli del parametro a.

Figura 13 : Evoluzione verso il Caos

La figura mostra in maniera grossolana l’evoluzione del sistema al

variare di a ed evidenzia l’intervallo dei valori nei quali il sistema si

comporta in maniera caotica. Superato questa finestra, il sistema diverge

all’infinito dato un qualsiasi stato iniziale.

Figura 14 : Passaggio dal ciclo instabile di periodo 2 a quello stabile di periodo 4

Figura 15 : Ciclo di periodo 8 e comportamento caotico