Università degli Studi di Siena
Dipartimento di Ingegneria
dell’Informazione e Scienze Matematiche
DI
Lecture Notes on
COMPLEX DYNAMIC SYSTEMS
MSc on Computer and Automation Engineering
MSc in Ingegneria Gestionale
MSc in Matematica
Chiara Mocenni
{[email protected]}
http://www.dii.unisi.it/∼mocenni/compdynsys-teach.html
A.A. 2013-2014
Indice
Introduzione
1
1 I sistemi lineari
2
2 Mappe e Punti Fissi
2.1 La mappa logistica . . . . . . . . . . . . . . . . . . . . . . . . .
2.2 Mappe di Poincaré . . . . . . . . . . . . . . . . . . . . . . . . .
2.3 Stabilità lineare di orbite periodiche . . . . . . . . . . . . . . . .
5
6
8
9
3 Biforcazioni di Sistemi Dinamici
3.1 Biforcazione Nodo-Sella . . . . . .
3.2 Biforcazione Transcritica . . . . . .
3.3 Biforcazione Pitchfork Supercritica
3.4 Biforcazione Pitchfork Subcritica .
3.5 Biforcazione di Hopf . . . . . . . .
4 Flussi sul cerchio e simulazioni
4.1 Pendolo semplice . . . . . . .
4.2 Pendolo sovrasmorzato . . . .
4.3 Beat phenomenon . . . . . . .
4.4 Lucciole . . . . . . . . . . . .
5 Il modello di Lorenz
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
10
11
12
13
14
15
.
.
.
.
18
18
21
22
24
32
6 Controllo del caos
42
6.1 Controllo di un processo caotico in continua evoluzione . . . . . 43
6.2 Targeting . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 44
6.3 Controllo del sistema di Lorenz . . . . . . . . . . . . . . . . . . 45
Bibliografia
54
i
Introduzione
Questo corso ha come obiettivo la presentazione delle basi metodologiche che
hanno portato alla definizione e all’analisi dei sistemi dinamici complessi. Tali metodologie trovano il loro fondamento nell’insieme delle teorie matematiche che consentono di modellizzare e studiare i sistemi dinamici non lineari,
come la teoria qualitativa delle equazioni differenziali, l’analisi della stabilià
delle soluzioni di equilibrio, la teoria delle biforcazioni e la teoria del caos
deterministico.
La ricerca di tali strumenti è stata stimolata dall’esigenza di rappresentare,
mediante modelli matematici, sistemi reali complessi che evolvono nel tempo,
come il moto dei pianeti, le oscillazioni di un pendolo, il flusso delle correnti
atmosferiche, lo scorrere più o meno regolare dell’acqua in un fiume, il numero di insetti che anno dopo anno popolano una certa regione, l’andamento
giornaliero dei prezzi delle azioni nei mercati finanziari e cosı̀ via.
Ciò che accomuna tali sistemi è la presenza di oscillazioni e fluttuazioni che
lasciano intuire la presenza di fenomeni stocastici e irregolarità dovute a fattori
esogeni. In realtà, si è visto che esistono una serie di sistemi dinamici descritti
da equazioni deterministiche (cioè senza componenti esogene stocastiche) non
lineari relativamente semplici che possono dar luogo a soluzioni estremamente
irregolari e imprevedibili. Tale dinamica prende il nome di caos deterministico.
L’apparente contraddizione contenuta nel termine caos deterministico nasce
dal fatto che i modelli matematici di tipo deterministico vengono in genere
associati all’idea di fenomeni regolari, prevedibili, che si ripetono nel tempo,
mentre il termine caotico viene riferito a situazioni caratterizzate da assenza
di regole e da imprevedibilità.
Grazie all’insieme degli strumenti che consentono di descrivere e analizzare
i sistemi dinamici complessi siamo in grado oggi di cogliere le caratteristiche
principali che si trovano alla base dei fenomeni complessi che osserviamo ogni
giorno nella natura e nella nostra vita.
1
Capitolo 1
I sistemi lineari
Si consideri il sistema dinamico descritto dalla seguente equazione differenziale
ordinaria:
ẋ = f (x),
x ∈ Rn ,
(1.1)
è la derivata temporale dello stato x. Una soluzione di (1.1) è una
dove ẋ = dx
dt
funzione x(x0 , t) a valori in Rn dipendente da t e condizione iniziale x(0) = x0 .
Possiamo interpretare la parte destra di (1.1) come un flusso φt : U → Rn ,
U ⊆ Rn , dove φt (x) = φ(x, t) è una funzione smooth, che soddisfa la (1.1),
definita per ogni x in U e t in qualche intervallo I = (a, b) ⊆ R.
Se il sistema è lineare e autonomo, cioè non dipende esplicitamente dal
tempo t, f (x) = Ax, dove A è una matrice a coefficienti costanti.
In questo caso la soluzione della (1.1) è
x(x0 , t) = etA x0 ,
(1.2)
dove etA è la matrice definita come segue:
etA = [I + tA +
tn
t2 2
A + · · · + An + . . . ].
2!
n!
(1.3)
Una soluzione generale x(t) della (1.1) può essere ottenuta comme combinazione di n soluzioni indipendenti {x1 (t), . . . , xn (t)}:
x(t) =
n
X
cj xj (t),
j=1
dove le n costanti incognite possono essere calcolate dalle condizioni iniziali.
Se A ha n autovettori linearmente indipendenti v j , j = 1, . . . , n, allora
possiamo utilizzare come base dello spazio delle soluzioni le funzioni
xj (t) = eλj t v j ,
2
3
dove λj è l’autovalore associato all’autovettore v j . Per una coppia di autovalori
complessi λ̄j = αj ±iβj con autovettori v R ±iv I , possiamo prendere come coppia
di soluzioni reali indipendenti ad essi associate
xj = eαj t (v R cos(βt) − v I sin(βt)),
xj+1 = eαj t (v R cos(βt) + v I sin(βt)).
Denotiamo con
X(t) = [x1 (t), . . . , xn (t)]
la matrice fondamentale delle soluzioni che ha le n soluzioni nelle colonne. Le
colonne xj (t), j = 1, . . . , n di X(t) formano una base dello spazio delle soluzioni
della (1.1), e inoltre
etA = X(t)X(0)−1 .
Nel caso in cui gli autovalori abbiano molteplicità maggiore di uno, è
possibile trovare una trasformazione invertibile che diagonalizza A o almeno permette di scriverla in forma normale di Jordan. L’equazione (1.1) diventa
allora
ẏ = Jy,
dove J = T −1 AT e x = T y e le colonne sono gli autovettori generalizzati di A.
L’esponenziale di A diventa allora
etA = T etJ T −1 ,
dove gli esponenziali sono calcolati per le tre matrici 2 x 2 in forma di Jordan:
λt
λ1 0
e 1
0
tA
A=
e =
;
0 λ2
0 eλ2 t
α −β
tA
αt cos(βt) −sin(βt)
A=
e =e
;
β α
sin(βt) cos(βt)
λ 0
1
0
A=
etA = eλt
.
1 λ
t 1
4
Consideriamo ora il caso in cui n = 2
a
A=
c
e
b
.
d
Il determinante della equazione caratteristica è
λ2 − τ λ + ∆ = 0,
dove
τ = traccia(A) = a + d;
∆ = det(A) = ad − bc.
Abbiamo
λ1,2 =
τ±
√
τ 2 − 4∆
,
2
∆ = λ1 λ2 ,
τ = λ1 + λ2 .
Questo permette di classificare la stabilità delle soluzioni per ogni combinazione
degli autovalori e di visualizzare questa classificazione nel piano (∆, τ ) come
segue:
• ∆ < 0: gli autovalori hanno segno opposto. Il punto di equilibrio x = 0
è un punto di sella, cioè il sistema è asintoticamente instabile;
• ∆ > 0: gli autovalori hanno lo stesso segno o sono complessi coniugati. Il
punto di equilibrio è un nodo se τ 2 −4∆ > 0 e una spirale se τ 2 −4∆ < 0.
La parabola τ 2 − 4∆ = 0 è il caso limite in cui si trovano le stelle e i nodi
degeneri. La stabilità dei nodi e delle spirali è determinata da τ : se τ < 0
tutti e due gli autovalori hanno parte reale negativa, quindi il punto
di equilibrio è asintoticamente stabile; se invece τ > 0 si hanno nodi
e spirali instabili. I centri (soluzioni marginalmente stabili) si trovano
nell’asse τ = 0;
• ∆ = 0: c’è almeno un autovalore nullo, ciò significa che l’origine non è un
punto di equilibrio isolato (stabilità marginale): esiste almeno una linea
di punti di equilibrio o addirittura un piano di punti di equilibrio (per
τ = 0).
Esempio: Studiare il sistema ẋ = x + y; ẏ = 4x − 2y, (x0 , y0 ) = (2, −3).
Capitolo 2
Mappe e Punti Fissi
Una mappa è descritta da una equazione del tipo
xn+1 = f (xn ),
(2.1)
dove f è una funzione continua con derivata continua, f : R → R. Vengono chiamate mappe le equazioni alle differenze, le iterazioni, e in generale le
equazioni in cui il tempo è discreto.
Definizione. x∗ si dice punto fisso per la mappa (2.1) se:
f (x∗ ) = x∗ .
(2.2)
Per la (2.2), se per qualche n accade che xn = x∗ , allora xn+1 = f (xn ) =
f (x∗ ) = x∗ , cioè l’orbita di (2.1) rimane uguale a x∗ per tutte le iterazioni
successive.
Per studiare la stabilità di x∗ si considera un’orbita vicina xn = x∗ + ηn
e si analizza il fatto che essa sia attratta o respinta da x∗ stessa, cioè si si
va a vedere se la perturbazione ηn cresce o si riduce. Sostituendo nella (2.1)
otteniamo:
x∗ + ηn+1 = xn+1 = f (x∗ + ηn ) = f (x∗ ) + f 0 (x∗ )ηn + O(ηn2 ).
(2.3)
Poichè f (x∗ ) = x∗ , l’equazione (6.19) diventa
ηn+1 = f 0 (x∗ )ηn + O(ηn2 ).
(2.4)
Se trascuriamo i termini in η 2 otteniamo una equazione lineare in η con
autovalore, detto anche in questo caso moltiplicatore, λ = f 0 (x∗ ).
La soluzione della mappa lineare si trova iterando a partire da un valore
iniziale η0 :
ηn = λn η0 .
(2.5)
Si distinguono tre casi:
5
2.1. LA MAPPA LOGISTICA
6
i) |λ| = |f 0 (x∗ )| < 1. In questo caso ηn → 0 per n → ∞ e il punto fisso è
linearmente stabile.
ii) |λ| = |f 0 (x∗ )| > 1. In questo caso ηn → ∞ per n → ∞ e il punto fisso è
linearmente instabile.
iii) |λ| = |f 0 (x∗ )| = 1. In questo caso il punto fisso è marginalmente stabile
e la sua stabilità è determinata dai termini non lineari.
2.1
La mappa logistica
Esempio. Determinare la stabilità dei punti fissi della mappa logistica:
xn+1 = f (xn ) = rxn (1 − xn ).
(2.6)
Assumiamo x ≥ 0. I punti fissi della mappa corrispondono alle intersezioni
del grafico della funzione f (xn ) con la bisettrice del piano (xn , xn+1 ), cioè
soluzioni tali che
xn , xn+1 .
Nel caso della mappa2.6 i punti fissi sono: x∗1 = 0 e x∗2 = 1 − 1/r. x∗1 è un punto
fisso per ogni valore di r, mentre x∗2 lo è solo per r ≥ 1.
Poiché la mappa è non lineare, la stabilità delle soluzioni di equilibrio (punti
fissi), può essere analizzata effettuando una linearizzazione intorno ai punti
fissi. Tale stabilità dipende dunque, localmente, dal moltiplicatore f 0 (x∗ ) =
r − 2rx∗ .
• x∗1 = 0. Poichè f 0 (x∗1 ) = r, l’origine è stabile per r < 1 ed instabile per
r > 1.
• x∗2 . f 0 (x∗2 ) = 2 − r, quindi x∗2 è stabile per −1 < (2 − r) < 1, cioè per
1 < r < 3, e instabile per r > 3.
• Per r > 3 entrambi i punti fissi sono instabili.
• Per r = 3 f 0 (x∗2 ) = −1, cioè r raggiunge un valore critico (per il quale
la stabilità della soluzione di equilibrio del sistema linearizzato è marginalmente stabile) superato il quale il punto x∗2 perde la sua stabilità e
le traiettorie del sistema cominciano ad oscillare indefinitamente tra due
valori.
2.1. LA MAPPA LOGISTICA
7
Questo tipo di perdita di stabilità prende il nome di biforcazione flip.
Quello che succede al sistema quando avviene una biforcazione flip è che un
ciclo di periodo T perde la sua stabilità dando luogo ad un ciclo di periodo
2T . Per capire questo fatto osserviamo che un ciclo di periodo 2T è un ciclo
di periodo T per la seconda iterata della mappa f : f 2 (x) = f (f (x)). Poichè in
questo caso f (x) è un polinomio quadratico, f 2 (x) è un polinomio di quarto
grado che piò avere fino a 4 intersezioni con la bisettrice. La biforcazione flip
può manifestarsi in ogni tipo di mappe.
Il ciclo di periodo 2 rimane stabile per 3 < r < 3.449. Oltre questo valore
anch’esso perde la propria stabilità dando luogo ad un nuovo ciclo di periodo
4T . Per comprendere questo fenomeno basta analizzare la quarta iterata della
mappa f (x) o la seconda iterata della mappa f 2 (x).
Questo meccanismo si ripete per valori sempre crescenti di r, fino a che per
r > 3.57 si osserva un comportamento aperiodico denominato caos deterministico.
E’ interessante osservare che la mappa logistica ha la seguente proprietà:
limn→∞
rn − rn−1
= 4.669,
rn+1 − rn
cioè esiste una costante universale δ = 4.669 chiamata costante di Feigenbaum che fissa il rapporto che lega valori consecutivi della costante di
biforcazione.
Il diagramma di biforcazione della mappa logistica è mostrato nella seguente figura:
Figura 2.1: Diagramma di biforcazione della mappa logistica.
2.2. MAPPE DI POINCARÉ
2.2
8
Mappe di Poincaré
Una classe molto importante di mappe è rappresentata dalle mappe di Poincaré. Sia S una superficie di dimensione n − 1 trasversale al flusso definito da
ẋ = f (x). Sia inoltre xk la k − esima intersezione del flusso con la superficie
S. Si definisce mappa di Poincaré la mappa seguente:
xk+1 = P (xk ),
dove P : S → S rappresenta la successione delle intersezioni successive del
flusso con S (vedi figura 2.2).
Figura 2.2: Mappa di Poincaré (source: Scholarpedia).
Per definizione, se x∗ è un punto fisso della mappa di Poincaré, allora
P (x∗ ) = x∗ . Ciò significa che la traiettoria del sistema è una curva chiusa.
Le mappe di Poincaré riducono il problema dell’esistenza di una curva chiusa per un sistema dinamico all’esistenza di un punto fisso per una mappa di
dimensione inferiore.
Si consideri ad esempio il sistema ṙ = r(1 − r2 ); θ̇ = 1. Se prendiamo S
uguale all’asse x positiva, il primo ritorno
R r1 sudr S avviene
R 2π dopo t = 2π. La mappa
di Poincaré r1 = P (r0 ) è definita da r0 r(1−r2 ) = 0 dt = 2π.
2.3. STABILITÀ LINEARE DI ORBITE PERIODICHE
2.3
9
Stabilità lineare di orbite periodiche
Sia ẋ = f (x) un sistema dinamico con un’orbita chiusa. Ci poniamo la domanda se questa orbita sia stabile, ovvero se il punto fisso x∗ della mappa di
Poincaré sia stabile. Si applica una perturbazione infinitesima v0 ad x∗ tale che
x∗ + v0 sia in S.
Dopo il primo ritorno si ha che:
x∗ + v1 = P (x∗ + v0 ) = P (x∗ ) + [DP (x∗ )] v0 + O(||v0 ||2 ),
dove DP (x∗ ) è la matrice della mappa di Poincaré linearizzata intorno a x∗ .
Allora la stabilità del punto fisso x∗ si esprime in base agli autovalori λj di
DP (x∗ ). L’orbita chiusa è stabile sse |λj | < 1, ∀j = 1, . . . , n − 1. Gli autovalori
λj prendono anche il nome di Moltiplicatori di Floquet.
Si consideri l’esempio precedente ṙ = r(1 − r2 ); θ̇ = 1. Sappiamo che esiste
un ciclo limite r∗ = 1. Sia r = 1 + η con η infinitesimo. Allora
ṙ = η̇ = (1 + η)(1 − (1 + η)2 ).
Se trascuriamo i termini non lineari, otteniamo
η̇ = −2η,
cioè
η(t) = η0 e−2t .
Dopo 2π, η1 = e−4π η0 e e−4π è il moltiplicatore di Floquet. Poiché |e−4π | < 1
si deduce che il ciclo limite r∗ = 1 è stabile.
Capitolo 3
Biforcazioni di Sistemi Dinamici
Un punto di equilibrio di un sistema dinamico di ordine n è detto iperbolico
se tutti gli autovalori del sistema linearizzato intorno al punto di equilibrio si
trovano fuori dall’asse immaginario.
L’importante teorema di Hartman-Grobman stabilisce che il phase portrait
(cioè il diagramma di tutte le traiettorie del sistema dinamico nello spazio delle fasi) vicino ad un punto iperbolico è topologicamente equivalente al phase
portrait del sistema linearizzato. Topologicamente equivalente significa che esiste un omeomorfismo (cioè una deformazione continua con inversa continua)
che trasforma un phase portrait locale nell’altro, mappa traiettorie in traiettorie e mantiene la direzione temporale.
In particolare, le proprietà di stabilità del punto di equilibrio sono completamente caratterizzate da quella della soluzione del sistema linearizzato.
Un phase portrait viene detto strutturalmente stabile se la sua topologia non
può essere modificata da una perturbazione arbitrariamente piccola del campo
vettoriale.
Se analizziamo le proprietà di un sistema dinamico in funzione dei parametri
che compaiono nelle equazioni, può accadere che vi siano dei cambiamenti nella
struttura qualitativa delle soluzioni per qualche valore dei parametri. Questi
cambiamenti vengono chiamati biforcazioni e i parametri che danno loro luogo
prendono il nome di valori di biforcazione.
Definizione 1. Dato il sistema dinamico
ẋ = fµ (x),
x ∈ Rn ,
µ ∈ Rk ,
(3.1)
un valore µ0 di (3.1) per cui il flusso definito dal sistema (3.1) non è strutturalmente stabile è un valore di biforcazione di µ.
Definizione 2. Il luogo dei punti nel piano (µ, x) degli insiemi invarianti di
(3.1) si chiama diagramma di biforcazione.
10
3.1. BIFORCAZIONE NODO-SELLA
11
Di seguito vengono elencati i tipi più comuni di biforcazioni locali in termini
della loro forma normale. La forma normale di una biforcazione è una rappresentazione prototipale di tutte le biforcazioni di quel tipo. Per chiarezza
esaminiamo sistemi di ordine 1. Nelle equazioni che seguono r è un parametro
che può essere positivo, negativo o nullo.
3.1
Biforcazione Nodo-Sella
ẋ = r + x2 .
(3.2)
20
15
10
5
0
−5
−10
−15
−6
−4
−2
0
2
4
Figura 3.1: Biforcazione Nodo-Sella. r0 = 0.
6
3.2. BIFORCAZIONE TRANSCRITICA
12
Figura 3.2: Diagramma di biforcazione: biforcazione nodo-sella.
3.2
Biforcazione Transcritica
ẋ = rx − x2 .
(3.3)
2
1.5
1
0.5
0
−0.5
−1
−1.5
−2
−2.5
−3
−3
−2
−1
0
1
2
3
Figura 3.3: Biforcazione Transcritica. r0 = 0.
3.3. BIFORCAZIONE PITCHFORK SUPERCRITICA
13
Figura 3.4: Diagramma di biforcazione: biforcazione transcritica.
3.3
Biforcazione Pitchfork Supercritica
ẋ = rx − x3 .
3.4
Biforcazione Pitchfork Subcritica
ẋ = rx + x3 .
3.5
(3.4)
(3.5)
Biforcazione di Hopf
In sistemi di ordine almeno 2 si incontra anche un altro tipo di biforcazione: la Biforcazione di Hopf. Supponiamo che il sistema linearizzato abbia
una coppia di autovalori complessi coniugati per il valore critico µ0 = 0 del
parametro µ:
λ = µ + iω, λ̄ = µ − iω
3.5. BIFORCAZIONE DI HOPF
14
3
2
1
0
−1
−2
−3
−2
−1.5
−1
−0.5
0
0.5
1
1.5
2
Figura 3.5: Biforcazione Pitchfork Supercritica. r0 = 0, α = −1, ω = 1, a =
−1, b = 0.
Figura 3.6: Diagramma di biforcazione: biforcazione pitchfork supercritica.
La forma normale della biforcazione di Hopf è la seguente:
ẋ = µx − ωy + (−x − by)(x2 + y 2 ) + O(|x|5 , |y|5 ),
ẏ = ωx + µy + (bx − y)(x2 + y 2 ) + O(|x|5 , |y|5 ).
(3.6)
3.5. BIFORCAZIONE DI HOPF
15
3
2
1
0
−1
−2
−3
−2
−1.5
−1
−0.5
0
0.5
1
1.5
2
Figura 3.7: Biforcazione Pitchfork Subcritica. r0 = 0.
Figura 3.8: Diagramma di biforcazione: biforcazione pitchfork subcritica.
L’equazione (3.6) in coordinate polari (x = rcos(θ), y = rsin(θ)) diventa:
ṙ = µr − r3 + O(|r|5 ),
θ̇ = ω + br2 + O(|r|4 ).
(3.7)
3.5. BIFORCAZIONE DI HOPF
16
x ’ = mu x − omega y + ( − x − b y) (x2 + y2)
y ’ = omega x + mu y + (b x − y) (x2 + y2)
mu = − 1 omega = 1
b=1
3
2
y
1
0
−1
−2
−3
−3
−2
−1
0
x
1
2
3
Figura 3.9: Biforcazione di Hopf Supercritica. µ < 0.
x ’ = mu x − omega y + ( − x − b y) (x2 + y2)
y ’ = omega x + mu y + (b x − y) (x2 + y2)
mu = 0
omega = 1
b=1
3
2
y
1
0
−1
−2
−3
−3
−2
−1
0
x
1
2
3
Figura 3.10: Biforcazione di Hopf Supercritica. µ = 0.
3.5. BIFORCAZIONE DI HOPF
17
x ’ = mu x − omega y + ( − x − b y) (x2 + y2)
y ’ = omega x + mu y + (b x − y) (x2 + y2)
mu = 1
omega = 1
b=1
3
2
y
1
0
−1
−2
−3
−3
−2
−1
0
x
1
2
3
Figura 3.11: Biforcazione di Hopf Supercritica. µ > 0.
Figura 3.12: Diagramma di biforcazione: biforcazione di Hopf supercritica.
Figura 3.13: Diagramma di biforcazione: biforcazione di Hopf subcritica.
Capitolo 4
Flussi sul cerchio e simulazioni
In questo capitolo viene descritta la la dinamica di sistemi oscillanti, quali il
pendolo semplice ed il pendolo sovrasmorzato. Si introduce inoltre il concetto
di sincronizzazione di sistemi dinamici.
4.1
Pendolo semplice
Si consideri il sistema pendolo governato dalle seguenti equazioni dinamiche:
θ̇ = ω
ω̇ = −sin(θ) − Dω
Con θ si rappresenta la posizione angolare del pendolo e con ω si intende la
velocità angolare dello stesso (in realtà è la velocità angolare del centro di massa
del sistema fisico pendolo). Il parametro D è la costante di proporzionalità
relativa all’attrito viscoso del sistema, infatti la componente Dω dipende dalla
velocità angolare del pendolo. Si considerino due casi distinti da analizzare nel
piano delle fasi. Il primo in cui D = 0 e il secondo in cui D 6= 0 (ma positivo).
Il piano delle fasi è caratterizzato da un campo vettoriale periodico di periodo
2π. Questo significa che si possono considerare due punti di equilibrio adiacenti
per descrivere il comportamento del piano delle fasi.
Caso 1: D = 0. Si consideri solo una finestra del piano delle fasi compresa
fra i punti (θ,ω)=(−π,0) e (2π,0), riportata in Figura 4.1. In questo spazio vi
sono quattro punti di equilibrio. Il primo, di coordinate x0 =(θ,ω)=(−π,0), è un
punto di equilibrio instabile perchè non è nè un attrattore nè Lyapunov stabile
(Re(λ) > 0). Il secondo punto, di coordinate x1 =(θ,ω)=(0,0), è un punto di
equilirio marginalmente stabile (Re(λ) = 0), in quanto per condizioni iniziali
18
4.1. PENDOLO SEMPLICE
19
nell’intorno di x1 , le traiettorie tendono a formare un centro in x1 . Quindi il
punto x1 è stabile ma non attrattivo (oscillatore armonico). Il terzo punto, di
coordinate x2 =(θ,ω)=(π,0), è un punto di equilibrio instabile (è un punto di
sella) perchè per una condizione iniziale nell’intorno della posizione di equilibrio
x2 , la traiettoria associata diverge dal punto di equilibrio. Dato che il sistema
non è smorzato (D = 0), una condizione iniziale nell’intorno della posizione di
equilibrio instabile fa si che la traiettoria ad essa associata, per t → ∞, tenda
a raggiungere la posizione di equilibrio instabile adiacente, oppure a portare il
sistema in oscillazione nell’intorno della posizione di equilibrio stabile pi vicina
(dipende dalla velocità iniziale). Fisicamente, considerando un filo inestensibile
a cui è collegata una massa, libero di ruotare attorno ad un fulcro, portando
il pendolo in posizione di equilibrio instabile e perturbando la sua posizione
iniziale (velocità non nulla), esso tenderà ad allontanarsi dal punto di equilibrio
instabile e continuerà ad oscillare per t → ∞. Il quarto punto, di coordinate
x3 =(2π,0) è un punto di equilibrio con caratteristiche uguali al secondo punto.
θ ’ =ω
ω ’ = − sin(θ) − D ω
4
D=0
t tendente a −inf
3
2
ω
1
0
−1
P.to instabile (sella)
−2
−3
P.to L.−stabile (centro)
Centro
−4
−3
−2
−1
0
1
2
3
4
5
6
θ
Figura 4.1: Piano delle fasi del pendolo: D = 0.
Caso 2: D 6= 0. Si sceglie D = 0.25 per la simulazione. In riferimento
4.1. PENDOLO SEMPLICE
20
alla Figura 4.1, il piano delle fasi risultante è diverso dal precedente per ciò
che riguarda i punti di equilibrio stabili. Infatti considerando una condizione
iniziale nell’intorno di un punto di equilibrio instabile, la traiettoria associata
converge verso il punto di equilibrio stabile adiacente formando una spirale
(λ ∈ C). In altri termini, la velocità del pendolo diminuisce gradualmente
per t → ∞ e la posizione angolare si assesta nell’intorno della posizione di
equilibrio stabile, al contrario di quanto accade per D = 0, caso in cui il
pendolo oscilla sempre con la stessa ampiezza (oscillatore armonico). I punti
di equilibrio stabili sono anche attrattori perchè le traiettorie con condizione
iniziale nell’intorno dei punti di equilibrio stabili convergono ad essi per t →
∞. Aumentando il parametro D si nota che diminuisce il tempo con cui una
traiettoria giunge nel punto di equilibrio stabile, partendo da un intorno di
esso, per t → ∞ e con movimento a spirale.
θ ’ =ω
ω ’ = − sin(θ) − D ω
4
D = 0.25
t tendente a −inf
P.to instabile (sella)
Spirale
3
2
ω
1
0
−1
−2
−3
−4
P.to L.−stabile e attrattore (spirale)
−3
−2
−1
0
1
2
3
4
5
6
θ
Figura 4.2: Piano delle fasi del pendolo: D = 0.25.
4.2. PENDOLO SOVRASMORZATO
4.2
21
Pendolo sovrasmorzato
In questa sezione si analizza un sistema oscillante al variare dei parametri
caratteristici, mettendo in luce la perdita di uniformità dinamica del sistema
stesso, cioè il moto del sistema non è più simmetrico.
Descrizione del modello
Si consideri il modello del pendolo fisico descritto dalla seguente equazione
differenziale del secondo ordine:
mL2 θ̈ + bθ̇ + mgLsin(θ) = Γ
(4.1)
in cui ogni parametro è positivo. In particolare m è la massa del pendolo, L
è la lunghezza dello stesso, b è il coefficiente di attrito viscoso del sistema e Γ è
una coppia costante positiva. Per inserire il modello nell’interfaccia di pplane7
si può riscrivere l’equazione differenziale del secondo ordine come sistema di
equazioni del primo ordine come segue:
θ̇ = ω
ω̇ = −
b
F
gsin(θ)
−
ω
+
L
mL2
mL2
A questo punto si considera il caso in cui il pendolo sia sovrasmorzato
imponendo b grande e il modello (4.1) diventa il seguente:
Γ
b
θ̇ =
sin(θ)
mgL
mgL
dove
τ=
mgL
Γ
t; γ =
b
mgL
Quindi si ottiene:
θ0 = γ − sin(θ)
dove
θ0 =
∂θ
∂τ
Simulazione
Gli angoli sono misurati in senso antiorario. La posizione di equilibrio instabile
vale π e quella stabile vale 0. Si analizzano i piani delle fasi risultanti al variare
di γ che rappresenta il rapporto fra la coppia esterna applicata e la massima
4.3. BEAT PHENOMENON
22
coppia gravitazionale.
Caso 1: γ → 1+ . Il piano delle fasi è mostrato in Figura 4.3.a. Sono visibili
anche le nullclines, le quali non si intersecano mai, dati i parametri associati
al modello (4.1), ma sono molto vicine ad ogni periodo della sinusoide in arancione (θ̇ = 0). Non esistono punti di equilibrio (F non viene mai annichilita),
pertanto per ogni condizione iniziale il pendolo ruota per t → ∞. Il moto del
pendolo è più rapido in discesa, partendo da una condizione iniziale pari a π
per raggiungere 2π, piuttosto che in salita, partendo da 0 verso π (ω >> 0).
Infatti in discesa agiscono sul pendolo sia la coppia F che la coppia gravitazionale, mentre in salita l’energia fornita dalla coppia F viene in parte persa con
l’azione opponente della coppia gravitazionale e dell’attrito viscoso (supposto
elevato). La Figura 4.3.b mostra una soluzione particolare del piano delle fasi
in funzione del tempo. In particolare con b = 5 e γ → 1+ .
Caso 2: γ > 1. Più è elevato il rapporto γ e più le nullclines si distanziano.
Non esistono punti di equilibrio. In questo caso il tempo di salita e il tempo di
discesa sono pressochè uguali, in quanto la coppia F è cosı̀ elevata da rendere
quasi trascurabili i termini relativi alla forza peso e all’attrito viscoso. Caso
3: γ = 1. Appare un punto di equilibrio in posizione θ = π/2. In Figura 4.4 è
rappresentato il piano delle fasi con γ = 1. Le nullclines si intersecano ad ogni
periodo della sinusoide ed è evidenziato in rosso il punto di equilibrio (nodal
skin) periodico.
Caso 4: γ < 1. In questo caso i punti di equilibrio sono due e sono tanto più
distanti quanto γ → 0. Il primo è uno spiral skin e il secondo è una sella. Il
piano delle fasi è mostrato in Figura 4.5.
Conclusioni. Si nota come al variare del rapporto γ si possa osservare un
sistema oscillante non uniforme. Partendo da γ = 0 (caso non contemplato
in quanto il sistema si riduce ad un oscillatore con due punti di equilibrio di
cui il minore è stabile, vedi il Capitolo 1), per γ → 1.5 si osserva come i due
punti di equilibrio si concentrino in uno solo (per γ = 1) per poi svanire completamente. L’analisi qui presentata può essere fatta anche considerando un
campo vettoriale sul cerchio, considerando un sistema oscillante non uniforme
semplice (di primo grado) e considerando il flusso sul cerchio, piuttosto che su
di una retta.
4.3
Beat phenomenon
In questo esempio si mostra un primo esempio di sincronizzazione fra sistemi.
4.3. BEAT PHENOMENON
23
Descrizione del modello
Nei sistemi oscillanti con due variabili indipendenti è interessante capire quando le due variabili assumono lo stesso valore. Si consideri l’esempio riportato
nel libro di testo, in cui vi sono due corridori che partono insieme su un circuito circolare. Il primo impiega T1 secondi per compiere un giro, mentre il
secondo impiega T2 > T1 secondi (è più lento). Bisogna capire quando il primo
corridore sorpassa il secondo. Dato che la velocità dei due corridori è costante,
i sistemi che li rappresentano sono descritti dalle equazioni:
θ̇i = ωi
dove i = 1, 2 sono i due corridori e ωi = 2π/Ti è la pulsazione naturale dei due
sistemi. Definendo la differenza di fase φ = θ1 − θ2 , si vuole capire quando essa
vale 2π, cioè quando il primo corridore, al tempo t, ha compiuto un giro in più
rispetto al secondo corridore. Quindi
φ̇ = θ̇1 − θ̇1 = ω1 − ω2 = 2π
In questo caso θ1 è in fase con θ2 ogni
T =
T1 T2
2π
=
ω1 − ω2
T1 − T2
Simulazione
Per la simulazione di questo sistema si utilizza un codice Matlab creato appositamente. Il sistema viene riscritto utilizzando il sistema di riferimento in
coordinate polari:
x
cos(ωt)
=
y
sin(ωt)
Si associa ad ogni corridore un sistema in coordinate polari. Per la simulazione
si sceglie la seguente configurazione:
• ω1 = 2;
• ω2 = 1;
Ciò significa che il primo corridore doppierà il secondo ogni
Tlap =
2π
= 2π
ω1 − ω2
4.4. LUCCIOLE
24
La condizione iniziale è (x1,2 , y1,2 ) = (0, 0). Si considera uno spazio temporale
di quattro periodi e si tracciano le curve φx = x1 − x2 e φy = y1 − y2 dei due
sistemi. Come mostrato in Figura 4.6, il primo corridore doppierà il secondo
ogni volta che (φx , φy ) = (0, 0), cioè ogni volta che i due grafici si intersecano
contemporaneamente sulla retta y = 0.
4.4
Lucciole
In questa sezione si analizza un secondo esempio di sincronizzazione tra sistemi.
Descrizione del modello
In correlazione con l’esempio precedente, vi è il fenomeno naturale della luce
periodica emessa dalle lucciole. La differenza sostanziale è che in questo caso, oltre al fatto che la luce emessa dalla lucciola è fortemente collegata alla
presenza di uno stimolo esterno, la differenza di fase tra l’emissione di luce e
lo stimolo non aumenta linearmente. Di seguito si fa una traduzione di alcune
parti del libro di testo per una migliore comprensione .
Modello. Supponendo che θ(t) sia la fase con cui la lucciola si illumina, dove
θ = 0 corrisponde al momento in cui viene emessa la luce, e che in assenza di
stimoli esterni la lucciola emette luce a frequenza ω, si può scrivere:
Θ̇ = Ω
dove Θ è un generico stimolo periodico. La risposta della lucciola può essere
modellata come segue: se lo stimolo è avanti nel ciclo, allora si assume che
la lucciola aumenti la frequenza di emissione per sincronizzarsi; se la lucciola
emette luce a una frequenza maggiore dello stimolo, allora tenterà di diminuire
la sua frequenza. Un modello semplice è il seguente:
θ̇ = ω + Asin(Θ − θ)
in cui A > 0. Per esempio, se Θ è avanti rispetto a θ (0 < Θ − θ < π) allora la
lucciola aumenta la frequenza di emissione di luce (θ̇ > ω). Il parametro A è
una misura dell’abilità della lucciola di modificare la sua frequenza istantanea
(resetting strength).
Analisi. Si consideri la differenza di fase φ:
φ̇ = Θ̇ − θ̇ = Ω − ω − Asin(φ)
4.4. LUCCIOLE
25
In questa forma è un oscillatore non uniforme. Come per il pendolo, si introducono due parametri che rendono il sistema adimensionale:
τ = At; µ =
Ω−ω
A
Allora
φ0 = µ − sin(φ)
dove φ0 = ∂φ/∂τ .
Simulazione
A questo punto si simula il modello con pplane7 al variare del parametro µ.
Caso 1: µ = 0. Come mostra la Figura 4.7, le traiettorie tendono a zero per
t → ∞ (in verde). Infatti φ = 0 è un punto di equilibrio in cui si concentrano
tutte le traiettorie con condizione iniziale nell’intorno di zero. Invece φ = π/2
è un punto di equilibrio instabile, in quanto le traiettorie tendono a divergere
per t → ∞ (in rosso).
Caso 2: 0 < µ < 1. I punti di equilibrio si avvicinano, ma mantengono le
stesse caratteristiche di stabilità, come mostra la Figura 4.8. Le nullclines del
sistema lucciola-stimolo si allontanano fra loro. In realtà il fatto che il punto
di equilibrio stabile φ∗ > 0 significa che la lucciola e lo stimolo esterno hanno
la stessa frequenza ma sono sfasati (phase-locked ).
Caso 3: µ = 1. In Figura 4.9 si nota come le traiettorie convergano all’unico
punto di equilibrio se hanno condizione iniziale minore di esso e divergano se
hanno condizione iniziale maggiore dello stesso. Il punto di equilibrio è mezzostabile.
Caso 4: µ > 1. In questo caso le traiettorie divergono per t → ∞ per qualunque condizione iniziale (Figura4.10. Non esistono punti di equilibrio e le
nullclines del sistema lucciola-stimolo non si intersecano mai. L’emissione di
luce della lucciola non sarà mai in fase con l’emissione di luce esterna. Questo
è dovuto al fatto che lo stimolo esterno ha una frequenza elevata e il sistema
lucciola non riesce ad agganciarsi in fase ad esso.
4.4. LUCCIOLE
26
θ ’ =ω
ω ’ = − g/L sin(θ) − b/(m L2) ω + F/(m L2)
b=5
m=1
F = 10
g = 9.81
L=1
2
3
4
4
3
2
ω
1
0
−1
−2
−3
−4
−5
−4
−3
−2
−1
0
θ
1
5
(a) Piano delle fasi
θ ’ =ω
ω ’ = − g/L sin(θ) − b/(m L2) ω + F/(m L2)
b=5
m=1
F = 10
g = 9.81
L=1
18
θ
16
ω
14
Discesa (rapida)
12
θ and ω
10
Salita (lenta)
8
6
X: 16.15
Y: 6.337
X: 15.07
Y: 3.177
4
2
0
0
2
4
6
8
10
t
12
14
16
18
(b) Soluzione particolare in funzione del tempo
Figura 4.3: Piano delle fasi con γ → 1+
20
4.4. LUCCIOLE
27
θ ’ =ω
ω ’ = − g/L sin(θ) − b/(m L2) ω + F/(m L2)
b=5
m=1
F = 9.81
g = 9.81
L=1
4
3
2
ω
1
0
−1
−2
−3
−4
−10
−8
−6
−4
−2
0
θ
2
4
6
Figura 4.4: Piano delle fasi con γ = 1
8
10
4.4. LUCCIOLE
28
θ ’ =ω
ω ’ = − g/L sin(θ) − b/(m L2) ω + F/(m L2)
b=5
m = 1 g = 9.81
F=5 L=1
2
3
4
3
2
ω
1
0
−1
−2
−3
−4
−5
−4
−3
−2
−1
0
θ
1
Figura 4.5: Piano delle fasi con γ < 1
4
5
4.4. LUCCIOLE
29
2
1.5
1
C.I.=(0,0)
φx , φy
0.5
(φ ,φ )=(0,0)
(φ ,φ )=(0,0)
x y
x y
0
X: 6.293
Y: 0.009814
−0.5
−1
−1.5
0
5
10
15
t
Figura 4.6: Beat phenomenon
φ ’ = O − o − A sin(φ)
O=1
A=1
o=1
4
Instabile
3
phi
2
1
0
Stabile
−1
−2
−2
−1
0
1
2
t
Figura 4.7: Lucciole: µ = 0
3
4
4.4. LUCCIOLE
30
φ ’ = O − o − A sin(φ)
O=1
A=1
o = 0.5
4
Instabile
3
2
Stabile
phi
1
0
−1
−2
−3
−4
−2
0
2
4
t
6
8
10
Figura 4.8: Lucciole: 0 < µ < 1
φ ’ = O − o − A sin(φ)
O=2
A=1
o=1
4
Instabile
3
2
phi
1
Stabile
0
−1
−2
−3
−4
−2
0
2
4
t
6
Figura 4.9: Lucciole: µ = 0
8
10
4.4. LUCCIOLE
31
φ ’ = O − o − A sin(φ)
O=3
A=1
o=1
4
3
2
phi
1
0
−1
−2
−3
−4
−2
0
2
4
t
6
Figura 4.10: Lucciole: µ > 1
8
10
Capitolo 5
Il modello di Lorenz
In un articolo del 1963, Deterministic Nonperiodic Flow, Lorenz cita le equazioni di Saltzman (1962), derivate dall’analisi di un sistema composto da un
fluido di profondità (altezza) uniforme, con differenti temperature tra la parte
superiore e la parte inferiore, in particolare con una variazione di temperatura lineare. Successivamente Rayleigh scoprı̀ un punto critico per il quale le
equazioni di Saltzman mostrano un moto convettivo, basato sul numero di Rayleigh. In questo capitolo si analizza il sistema di Lorenz a partire dal modello
matematico.
Modello matematico
Lorenz definisce un set di tre equazioni dinamiche. In particolare x è proporzionale all’intensità del moto convettivo, y è proporzionale alla differenza di
temperatura tra le correnti ascensionali e discensionali e z è proporzionale alla distorsione del profilo verticale di temperatura (dalla linearità). Il sistema
tridimensionale è il seguente:
ẋ = −σx + σy
ẏ = −xz + rx − y
ż = xy − bz
dove, in questo caso, σ = 10 (numero di Prandtl ), r è il numero di Rayleigh
e b = 8/3 è la largezza approssimativa della regione di lavoro. Di seguito si
analizzano alcune peculiarità del sistema di Lorenz, con particolare attenzione
alle biforcazioni del sistema al variare del numero di Rayleigh 1 .
1
Il modello della convezione del fluido è realistico per valori di r compresi fra 0 e 1. Per
32
33
Analisi con matcont per 0 ≤ r < 30
In questa sezione si analizza empiricamente il sistema di Lorenz al variare del
parametro r, con σ e b fissati.
Caso 1: 0 ≤ r < 1. Il punto di equilibrio stabile, di coordinate (x, y, z) =
(0, 0, 0), è anche un attrattore globale. In Figura 5.1 si riportano tre simulazioni per r = 0; 0.5; 1 − ε. Per comodità viene rappresentata la variabile x(t)
in funzione del tempo.
10
8
6
4
x
2
0
−2
0 < r <1
−4
−6
−8
−10
0
0.5
1
1.5
2
2.5
t
3
3.5
4
4.5
5
Figura 5.1: Caso 1: 0 ≤ r < 1
Caso 2: 1 < r < 24. Per r = 1 si verifica la prima biforcazione del sistema.
In Figura 5.2 si può osservare come le traiettorie, per t → ∞, generate a
partire da condizioni iniziali opposte in segno, si allontanino definitivamente
dalla precedente soluzione di equilibrio stabile ((0, 0, 0)) e si distanzino fra loro
al variare di r. L’incremento del parametro r comporta la perdita di stabilità
della soluzione di equilibrio (0, 0, 0) e la generazione di due nuove soluzioni di
equilibrio stabili e simmetriche nel diagramma di biforcazione. Risulta chiaro
che la biforcazione è di tipo pitchfork. Nel caso riportato in Figura 5.3, in cui
valori maggiori di r si modellano sistemi come spiking irregolare nei lasers, convezione in
regioni toroidali, etc.
34
r = 22, le oscillazioni del sistema sono molto sostenute, ma comunque smorzate
(spirale). Inoltre, osservando il transitorio di x(t), si nota la presenza di uno
spike. Di seguito si comprende come mai vi sia questo fenomeno.
10
8
traiettoria sovrasmorzata
6
4
2
x
r incrementa
0
−2
−4
r è circa 1 e le traiettorie sono stabili verso (0,0,0)
−6
−8
−10
0
0.5
1
1.5
2
2.5
t
3
3.5
4
4.5
5
Figura 5.2: Caso 2: 1 < r < 24
Caso 3: 24 < r < 30. Per r = 24.7 si verifica la seconda biforcazione del
sistema. In realtà, per la simmetria della biforcazione pitchfork, le biforcazioni
sono due e sono simmetriche rispetto all’asse x = 0. Si nota che i due autovalori complessi e coniugati λ2 , λ3 hanno parte reale strettamente minore di
zero per r = 24.7 − δ e assumono parte reale strettamente maggiore di zero
per r = 24.7 + δ. In questo caso l’equilibrio è instabile e la biforcazione è di
Hopf subcritica. Le traiettorie tridimensionali, come enuncia la teoria classica,
dovrebbero allontanarsi definitivamente dagli autovettori instabili. In realtà
nessuna traiettoria diverge, ma si nota che esse entrano in due regioni diverse
senza un criterio intuitivo, con moto a spirale, senza mai stabilizzarsi nell’origine della spirale stessa. In riferimento alla Figura 5.5, si osserva il moto a
spirale caotica in tridimensione. Il segnale è assolutamente aperiodico, rendendo impredicibile il comportamento del sistema per valori futuri del tempo,
rispetto ad un tempo t0 generico. Da un sistema deterministico e totalmente
predicibile si passa ad un sistema caotico impredicibile e, come mostra la Figura 5.4, molto sensibile alle condizioni iniziali, per un r fissato. Infatti per una
piccolissima variazione (dell’ordine di 10−10 ) della condizione iniziale relativa
35
alla coordinata z si ottengono due risposte temporali differenti.
In Figura 5.6 è rappresentato il diagramma di biforcazione del modello di
Lorenz che riassume i tre casi analizzati.
Altri casi notevoli
Si consideri il caso in cui 144 < r < 166. In Figura 5.7 è rappresentata la risposta del sistema per r = 145. Si nota che le traiettorie tendono a stabilizzarsi su
un’orbita periodica (per r crescente), ma tornano ad essere aperiodiche dopo
un determinato intervallo di tempo. Questo fenomeno è chiamato caos intermittente e accade per una finestra di valori di r, tale per cui si verifica il period
doubling. Se r è dentro questo range, il caos intermittente può eventualmente
cessare, ma se r assume un valore prossimo all’ upper bound (r = 166.3), allora
il fenomeno si verifica per ogni valore di t. Infatti per r crescente da 145 a 166
si nota che il periodo in cui domina il caos intermittente tende ad aumentare
fino a dominare l’intera evoluzione temporale.
Nel caso in cui r = 212 si osserva un secondo fenomeno di period doubling,
ma con periodo raddoppiato rispetto alla finestra considerata in precedenza.
La risposta dinamica del sistema segue una traiettoria stabile, ma in questo
caso il caos non è di tipo intermittente, in quanto, per r crescente, la periodicità
della soluzione x(t) è di esattamente un periodo per ogni intervallo temporale.
In Figura 5.8 è rappresentata la realizzazione del sistema per r = 212.
Conclusioni
Dall’analisi proposta emergono alcune delle caratteristiche principali del sistema di Lorenz.
Aperiodicità. Per un determinato valore del numero di Rayleigh si osserva
il fenomeno dell’aperiodicità. E’ del tutto contro-intuitivo pensare che da un
sistema lineare si riesca a generare un segnale aperiodico, non predicibile e
apparentemente casuale (in un certo senso è anche non causale).
Sensitività. Il sistema risulta sensibile alle condizioni iniziali per quanto riguarda la risposta aperiodica. Un fattore di differenza molto piccolo genera
due segnali totalmente differenti dopo un tempo t. Dall’abstract dell’articolo
di Lorenz del 1963: [... soluzioni nonperiodiche sono instabili rispetto a piccole
perturbazioni, cosicchè condizioni iniziali di poco differenti possono evolvere
in stati considerevolmente diversi...].
Biforcazioni. Al variare di r vi sono tre biforcazioni di cui due sono simmetri-
36
che, ma le traiettorie, per qualsiasi valore di r, non divergono verso l’infinito.
Un sistema con questo comportamento si chiama dissipativo.
Simmetria. Il sistema è invariante se si inverte il segno delle variabili di stato
(tranne z). L’invarianza dell’asse z implica che tutte le traiettorie relative ad
esso restano sullo stesso e si avvicinano all’origine. Un fenomeno di simmetria
si osserva anche per le biforcazioni di Hopf, dal momento che si generano a
partire dai due rami stabili della biforcazione pitchfork.
37
20
15
10
0
−5
−10
−15
−20
0
5
10
15
t
20
25
(a) x(t)
40
35
30
25
20
z
x
5
15
10
5
0
0
5
10
15
20
x
−20
−15
−10
0
−5
5
y
(b) Plot tridimensionale
Figura 5.3: Caso 2: r=22
10
15
20
30
38
differenziazione traiettorie
si "gira" nella prima spirale
20
15
10
x
5
0
scambio
−5
−10
−15
si "gira" nella seconda spirale
−20
0
5
10
15
20
t
25
30
35
40
Figura 5.4: Risposte temporali per r = 27 e condizioni iniziali differenti per
10−10
40
35
30
z
25
20
15
20
10
15
5
10
5
0
20
15
10
5
0
0
−5
−10
−15
y
Figura 5.5: Moto caotico a spirale
−20
x
39
30
H
25
H
r
20
15
10
5
BP
0
−10
−8
−6
−4
−2
0
x
2
4
6
Figura 5.6: Biforcazioni al variare del numero di Rayleigh
8
10
40
40
30
20
x
10
0
−10
−20
−30
−40
0
2
4
6
t
8
10
12
(a) x(t)
250
200
z
150
100
50
0
300
200
100
0
−100
150
50
100
0
y
(b) Plot tridimensionale
Figura 5.7: Risposta del sistema per r = 145
−50
−100
41
250
200
z
150
100
400
50
300
200
0
100
150
100
0
50
0
−50
−100
−100
x
y
(a) Plot tridimensionale
40
30
20
10
x
0
−10
−20
−30
−40
−50
0
2
4
6
8
10
t
12
14
(b) x(t)
Figura 5.8: Risposta dinamica per r = 212
16
18
20
Capitolo 6
Controllo del caos
Per completare la panoramica sul caos deterministico, si analizza il problema
del controllo del caos. Si nota che al giorno d’oggi si utilizzano tecniche di
controllo che derivano senza alterazioni dalla teoria classica del controllo.
Introduzione
Un sistema caotico è caratterizzato dalla dipendenza dalle condizioni iniziali.
Due traiettorie che partono da punti molto vicini nello spazio divergono, ad
esempio, esponenzialmente. Infatti, in un processo di identificazione di un sistema (black-box ), spesso non si riesce a distinguere esattamente se si identifica
un sistema affetto da rumore oppure un sistema caotico. Un quantificatore della divergenza delle due traiettorie è l’esponente di Lyapunov, in particolare il
suo massimo. Un secondo problema è la possibilità di sincronizzazione fra due
sistemi caotici identici. Recenti studi dimostrano che effettivamente può avvenire la sincronizzazione, con importanti applicazioni come telecomunicazioni
digitali non lineari (Pecora e Carrol, 1990).
Il problema messo in luce in questo lavoro è la possibilità di stabilizzare
un sistema caotico (di Lorenz) con un semplice regolatore PI (proporzionaleintegrale).
Teoria del controllo
Con la parola controllo si intende controllo in retroazione. Si consideri un sistema non lineare. Si hanno alcune variabili dipendenti dal tempo e si può pensare
di controllarle, ossia di modificarne lo stato e la sua evoluzione, conoscendo la
42
6.1. CONTROLLO DI UN PROCESSO CAOTICO IN CONTINUA
EVOLUZIONE
43
storia del sistema e il suo stato attuale. Il controllo può essere inteso in tre
modi diversi, elencati di seguito:
1. Dato un sistema caotico in costante evoluzione, aumentare la sua performance media nel tempo, dove performance è un termine generico per
definire una caratteristica particolare che si vuole ottenere dal sistema;
2. Dato un sistema caotico con una condizione iniziale associata, trascinare
rapidamente la traiettoria rapidamente verso lo stato desiderato;
3. Dato un processo caotico in costante evoluzione, prevenire che diverga
verso una regione non buona dello spazio delle soluzioni.
L’esistenza di un’orbita complessa è molto rilevante per ottenere il target
assegnato nel primo dei tre problemi. Per il secondo e il terzo target è più importante la sensibilità esponenziale alle condizioni iniziali. Per tutti i targets è
da notare che il controllo può essere ottenuto soltanto con piccoli controlli, cioè
con piccole variazioni dello stato attuale. L’interesse per il terzo target risiede
nella possibilità di controllare e prevenire eventi catastrofici che si verificano
nel momento in cui lo stato converge in una precisa zona dello spazio delle
soluzioni.
Ovviamente esiste un approccio diverso per ognuno dei tre targets. Nelle
sezioni seguenti si tratta del primo e del secondo tipo di approccio rispettivamente per il primo e per il secondo target.
6.1
Controllo di un processo caotico in continua evoluzione
Il punto chiave per questo problema è la presenza di un’ orbita complessa nel
caos. Con questo si intende che esistono molte orbite topologicamente distinte
all’interno di un set caotico invariante. Quindi vi sono tipicamente molte orbite
periodiche instabili (UPO - Unstable periodic orbit) nascoste in un attrattore
caotico. Infatti, per un attrattore iperbolico, il numero di UPOs con periodo
minore o uguale di P incrementa esponenzialmente con P ,
#(P ) ∼ exp(hT P )
(6.1)
dove hT è l’entropia topologica (Katok, 1980). Per risolvere il problema esiste
l’approccio seguente (Ott et al. 1990a,b):
1. Esaminare le dinamiche non controllare del sistema e determinare alcune
UPOs di periodo piccolo nascoste nell’attrattore;
6.2. TARGETING
44
2. Esaminare queste UPOs e trovare la caratteristica risultante dalla dinamica temporale di ogni UPO. Confrontare le caratteristiche delle UPOs
e scegliere la migliore.
3. Formulare un algoritmo di controllo che stabilizzi la UPO selezionata in
un intorno di essa.
In sostanza si cerca di far convergere una traiettoria verso un’orbita propria del sistema e mantenere questo comportamento per tutta la durata del
processo. Nel caso in cui una traiettoria tenti di divergere da quest’orbita, si
cerca di intervenire tempestivamente per riportarla sull’orbita. La tempestività è importante perchè il sistema può essere controllato in maniera semplice
solo in un intorno molto ristretto della UPO selezionata. Si pensi all’attrattore
di Lorenz. Sostanzialmente esiste un’orbita in cui le traiettorie convergono e
successivamente divergono e si spostano verso una seconda orbita, tentando
di convergere, ma, per la naturale evoluzione del sistema non controllato di
Lorenz, divergono anche da essa e si sposteranno verso la prima orbita, reiniziando il ciclo. Se si identificasse una regione Ω con l’approccio di Ott, si
potrebbe pensare di stabilizzare la traiettoria nell’intorno di una delle due orbite caotiche. Data la natura del sistema non si può pensare di stabilizzare la
traiettoria nell’origine, in quanto non si saprebbe come governare il sistema
verso un punto che non è più di equilibrio (nel caso caotico). A questo punto
diventerebbe un problema di inseguimento e non più di stabilizzazione, ma in
realtà un sistema caotico non ha ingressi esterni (è autonomo di ed ha caos
autonomo).
Per identificare una UPO e le sue caratteristiche si considerano delle mappe ricavate sezionando con un piano un’orbita dell’attrattore e analizzando le
serie temporali derivate. In questo modo si trasferisce il problema tridimensionale in un problema bidimensionale. Ovviamente la riduzione del problema ad
un caso più semplice porta ad avere un’analisi semplificata. Nella fase in cui si
analizzano le serie temporali vi è in gioco il fattore rumore. Come precedentemente illustrato è un fattore di non facile identificazione, soprattutto quando
non si conosce il sistema da cui proviene l’attrattore strano. Per stabilizzare
una UPO si utilizza la teoria del controllo lineare. Si considerano quindi approcci come la tecnica del pole placement oppure un regolatore generico PI,
come descritto nell’esempio di questa relazione.
6.2
Targeting
Per la seconda tipologia di controllo del caos si analizza il problema di direzionare l’orbita di un sistema caotico in una zona desiderata del piano delle fasi.
6.3. CONTROLLO DEL SISTEMA DI LORENZ
45
Alla base del targeting vi è l’idea che un sistema caotico è esponenzialmente
sensibile a piccole perturbazioni. Una piccola perturbazione ha un largo effetto
in un tempo relativamente piccolo. Se le perturbazioni sono rumori ignoti, allora il sistema diventa impredicibile. Ma se il rumore presente è relativamente
piccolo, si possono tranquillamente inserire piccole perturbazioni per controllare il sistema. Una possibile strategia di controllo è il metodo forward-backward
(Shinbrot et al., 1990), descritto di seguito. Si consideri una mappa bidimensionale con una direzione di contrazione (un esponente di Lyapunov negetivo)
e una direzione di espansione (un esponente di Lyapunov positivo). Considerando uno scalare di controllo p, si può perturbare il punto di partenza verso
un set di regioni nello spazio delle soluzioni, cioè si cerca di indirizzare il punto
di partenza verso la regione desiderata. L’azione di controllo è manifestata nel
sistema come perturbazione iniziale, cioè il controllo è attivo solo inizialmente
(per t = 0). A questo punto si utilizza la conoscenza a priori del sistema, quindi
la conoscenza del set di traiettorie che il punto potrebbe percorrere evolvendo in un tempo T. Simultaneamente si fa evolvere la regione a cui si punta
indietro nel tempo T . Incrementando T , l’immagine del segmento di controllo
iniziale (in avanti nel tempo) e l’immagine della regione puntata (indietro nel
tempo) si espandono esponenzialmente in lunghezza. Nel momento in cui la
preimmagine della regione puntata e l’immagine del segmento di intersecano,
si identifica un valore del segmento iniziale p. Idealmente la traiettoria del
punto iniziale tenderà a collidere con la regione puntata in un tempo 2T . In
un sistema reale, affetto da rumore, si commette un errore anche molto piccolo quando si imposta il segmento iniziale: l’errore aumenta esponenzialmente
nel tempo ed e è probabile che il controllo adottato non sia robusto. Questo
sistema funzione per una dimensione del sistema non troppo elevata. Per sistemi di dimensione alta esiste il metodo di Kostelich (Kostelich et al., 1993).
La prima applicazione del controllo formalizzato da Kostelich si ha nel 1985
da parte della NASA per il targeting di una sonda spaziale con la cometa di
Giacobini-Zinner.
6.3
Controllo del sistema di Lorenz
In questa sezione si analizza la possibilità di controllare alcuni punti instabili
del sistema di Lorenz e stabilizzarli con la teoria del controllo lineare, in particolare con un regolatore proporzionale e integrale, utilizzando l’approccio di
Ott (primo metodo). Si consideri il sistema seguente:
ẋ = Ax + g(x) + u
(6.2)
6.3. CONTROLLO DEL SISTEMA DI LORENZ
46
dove x ∈ Rn è il vettore di stato, u ∈ Rn è il comando da progettare, A ∈
Rn×n è una matrice costante e g(x) è una funzione non lineare continuamente
differenziabile. Considerando xs un punto di equilibrio instabile del sistema,
tale che
Axs + g(xs ) = 0
(6.3)
si può progettare un controllore u che stabilizzi la traiettoria verso xs . Il
controllore u è un PI della forma
Z t
(6.4)
u = λ[B(Kx + k (y − ys )dτ )]
0
1×n
dove B ∈ R
è un vettore costante di guadagni, K ∈ R1 timesn è il
vettore di guadagni di retroazione proporzionale allo stato, k ∈ R è il guadagno
integrale, y = Cx è l’uscita con una matrice costante C ∈ R1×n , ys = Cxs è
l’osservazione del punto di equilibrio target xs , e λ ∈ {0, 1} è una variabile
binaria di attivazione che vale 1 se il vettore di stato ad un certo tempo t è
dentro la regione di controllo Ωxs di xs . Quindi il controllore u è attivo se e solo
se la traiettoria del punto x(t) è dentro la regione ΩxS . Si riscrive il sistema
come segue:
Z t
(6.5)
ẋ = Ax + g(x) + λ B Kx + k (Cx − Cxs )dτ
0
Si consideri il sistema di Lorenz. Imponendo le condizioni per i punti di equilibrio (derivate prime uguali
identificano tre punti di equilibrio
dipen-p
p
p a zero), sip
denti dai parametri: (− β(ρ − 1), − β(ρ − 1), (ρ−1)), (0, 0, 0) e β(ρ − 1), β(ρ − 1), (ρ−
1)). Riscrivendo il sistema nella forma 6.5, segue che l’errore di inseguimento
è descritto dalla seguente formula:


 

0
0
0
ex
0
g(x) − g(xs ) =  −zs 0 −xs   ey  +  −ex ez  = g 0 (xs )e + ḡ(e) (6.6)
y s xs 0
ez
ex ey
La dinamica dell’ errore incrementale è descritto dalle seguenti equazioni di
stato:
ė1 = (Ã + λB̄ K̄)e1 + g̃(e1 )
(6.7)
dove
e1 =
e
q
(6.8)


0
 −ex ez 
ḡ(e)

g̃(e1 ) =
=
 ex ey 
0
0
(6.9)
6.3. CONTROLLO DEL SISTEMA DI LORENZ
47


−σ
σ
0
à = A + g 0 (xs ) =  ρ − zs −1 −xs 
ys
xs −β


−σ
σ
0
0
 ρ − zs −1 −xs 0 
Ā 0

à =
=
 ys
C 0
xs −β 0 
c1
c2
c3 0


b1

B
b2 

B̄ =
=
 b3 
0
0
K̄ = k1 k2 k3 k
(6.10)
(6.11)
(6.12)
(6.13)
Se la matrice à + B̄ K̄ è stabile in xs , allora si può affermare che il controllore
stabilizza l’evoluzione del sistema verso il target (la 6.7 è asintoticamente
stabip
p
le). Se la stessa matrice è stabile nel punto di equilibrio P1 =(− β(ρ − 1), − β(ρ − 1), (ρ−
1))), allora la 6.7 è asintoticamente stabile e il sistema di Lorenz può essere
stabilizzato verso il punto
P1 . Le stesse osservazioni valgono per
p
p di equilibrio
P2 = (0, 0, 0) e P3 = ( β(ρ − 1), β(ρ − 1), (ρ − 1)). Ad esempio si può pensare di stabilizzare il sistema di Lorenz per σ = 16, ρ = 45.6 e β = 4, per i
quali il sistema ha un comportamento caotico. Il sistema è stabilizzabile con
il controllore PI verso il punto di equilibrio P3 perchè la matrice à + B̄ K̄
è stabile in xs ha autovalori a parte reale minore di zero, quindi è stabile
(λ1 = −2.8, λ2 = −2.7, λ3 = −2.3, λ4 = −2.5). In questo caso le matrici del
sistema di controllo sono le seguenti
 
1

B= 0 
(6.14)
0
C = −1 1 0
(6.15)
K = 10.70 −19.18 12.91 , k = 0.12
(6.16)
In Figura 6.1 è riportato lo schema Simulink realizzato ad hoc per la simulazione. Le traiettorie delle tre variabili di stato tendono asintoticamente
al punto di equilibrio P3 = (13.35, 13.35, 44.60) pertanto l’obiettivo proposto è
stato raggiunto. Nelle Figure 6.2, 6.3 e 6.4 si mostrano le traiettorie temporali
rispettivamente per x(t), y(t) e z(t).
6.3. CONTROLLO DEL SISTEMA DI LORENZ
48
Figura 6.1: Modello Simulink del sistema di Lorenz controllato con un PI
110
100
90
80
x
70
60
50
40
30
20
10
0
50
100
150
200
250
t
300
Figura 6.2: x(t)
350
400
450
500
6.3. CONTROLLO DEL SISTEMA DI LORENZ
49
20
15
y
10
5
0
−5
−10
0
50
100
150
200
250
x
300
Figura 6.3: y(t)
350
400
450
500
6.3. CONTROLLO DEL SISTEMA DI LORENZ
50
60
55
z
50
45
40
35
30
0
50
100
150
200
250
t
300
Figura 6.4: z(t)
350
400
450
500
Linee guida per l’analisi dei
sistemi dinamici caotici
Nel seguito vengono presentate alcune linee guida per effettuare l’analisi dei
sistemi dinamici non lineari che mostrano comportamenti complessi, come ad
esempio le dinamiche caotiche.
1. Soluzioni di equilibrio. In primo luogo è opportuno individuare le
soluzioni di equilibrio del sistema dinamico. Queste soluzioni sono tali
che:
ẋ = 0, x ∈ Rn , x = (x1 , x2 , . . . , xn ).
(6.17)
In altre parole si tratta di risolvere il sistema:
ẋ1 = 0
ẋ2 = 0
·
·
ẋn = 0.
Ad esempio, nel caso del sistema di Lorenz, descritto dal sistema di
equazioni
ẋ = −σx + σy
ẏ = −xz + rx − y
ż = xy − bz,
si hanno tre soluzioni di equilibrio ottenute risolvendo il sistema
−σx + σy = 0
−xz + rx − y = 0
xy − bz = 0.
51
(6.18)
6.3. CONTROLLO DEL SISTEMA DI LORENZ
52
Le soluzioni sono:
√ √
√ √
x∗1 = (0, 0, 0); x∗2,3 = (± b r − 1, ± b r − 1, r − 1).
E’ evidente che le soluzioni x∗2,3 esistono solo per r ≥ 1.
2. Linearizzazione. Si studia il sistema linearizzato nell’intorno di ciascuno dei punti di equilibrio. In particolare, si calcolano gli autovalori della
matrice Jacobiano J nel punto di equilibrio in questione allo scopo di
studiare la stabilità della soluzione nulla del sistema linearizzato:
η̇ = f 0 (x∗ )η,
(6.19)
dove η = x − x∗ e J(x∗ ) = f 0 (x∗ ).
3. Stabilità degli equilibri in funzione dei parametri. Nel caso in cui
si abbiano alcuni autovalori a parte reale nulla (cioè tali che i corrispondenti punti di equilibrio sono non iperbolici) si devono effettuare alcune
simulazioni variando le condizioni iniziali ed il parametro critico (cioè
quello per cui si annulla la parte reale degli autovalori) alternativamente
allo scopo di comprendere la dinamica del sistema. Le simulazioni possono essere fatte utilizzando una delle funzioni ODE di MATLAB, WINPP,
MATCONT o qualunque software di integrazione.
4. Analisi di biforcazione. Si effettua un’analisi di biforcazione del sistema. A questo scopo si prende in considerazione una delle soluzioni
di equilibrio stabili trovate precedentemente. Si pongono i valori di tutti i parametri corrispondenti ad una delle condizioni di stabilità e si
eseguono delle simulazioni imponendo piccole variazioni del parametro
individuato come critico. Eventualmente si può utilizzare un software di
continuazione (ad esempio MATCONT).
Nel caso del sistema di Lorenz si evidenzia quanto segue:
- 0 ≤ r < 1. L’unica soluzione di equilibrio x∗1 è stabile asintoticamente.
- 1 < r < 24. La soluzione x∗1 perde la propria stabilità attraverso una
biforcazione pitchfork. Le due nuove soluzioni che nascono grazie alla
stessa biforcazione x∗2,3 risultano stabili per questi valori dei parametri.
- 24 < r < 30. Le due soluzioni x∗2,3 perdono stabilità attraverso una biforcazione di Hopf subcritica. Questa biforcazione non può essere facilmente
detectata in quanto il ciclo limite che esiste prime della biforcazione è
instabile.
- r > 24.74 nessuna delle soluzioni di equilibrio esistenti è stabile. Si
6.3. CONTROLLO DEL SISTEMA DI LORENZ
53
osservi che per questo valore di biforcazione non si determina nessuna
collisione di equilibri o cicli e tutte le soluzioni di equilibrio continuano
ad esistere, come confermato dall’analisi effettuata al punto [1].
Bibliografia
[1] Edward Ott Chaos in dynamical systems Cambridge University Press,
1993.
[2] G. Jiang, G. Chen, W.K. Tang Stabilizing unstable equilibrium points
of a class of chaotic systems using a state PI regulator (IEEE Trans. on
Circuits and Systems, VOL. 49, NO. 12, 2002)
[3] Steven H. Strogatz Nonlinear Dynamics and Chaos Westview, 1994.
54