Alma Mater Studiorum · Universita di Bologna

FACOLTA DI SCIENZE MATEMATICHE, FISICHE E NATURALICorso di Laurea in Scienze di Internet

STABILITA DEI SISTEMI LINEARIE

APPLICAZIONI ECONOMICHE

Tesi di Laurea in Matematica Generale

Relatore:Chiar.mo Prof.RITA FIORESI

Presentata da:MATTEO MALAGOLI

Sessione IIIAnno Accademico 2007/2008

Introduzione

La tesi nasce con l’intento di studiare e modellare un sistema economico competitivo, cioecercheremo di analizzare la competizione al variare del tempo attraverso un approccioanalitico e computazionale.

Il primo capitolo introdurra molte tipologie di equazioni differenziali, tra cui leequazioni non lineari. Per risolvere queste particolari equazioni procederemo verso unapproccio di tipo geometrico e alla linearizzazione del sistema.Partendo da una descrizione sintetica delle fasi, procederemo alla caratterizzazione ditutte le differenti tipologie di punti critici. In un primo momento l’analisi vertera susistemi lineari, per poi approdare ai sistemi non lineari, ed infine allo studio della stabilitadei loro punti di equilibrio attraverso il metodo indiretto di Lyapunov. Infine studieremoil modello Competizione tra specie, modello base di riferimento per lo studio del nostroproblema economico.

Il secondo capitolo riprendera i concetti esaminati chiarendo i quesiti risolti e quellinon ancora risolti. Si comprendera l’impossibilta di reperire informazioni riguardanti ibacini di attrazione, ossia quei domini di grande importanza per la comprensione delsistema.Per sopperire alla mancanza di informazioni introdurremo il metodo diretto di Lyapunov,e i relativi risultati.

2

Nel terzo capitolo si utilizzeranno gran parte delle competenze apprese nei dueprecedenti applicate al Fenomeno Empirico. Si introdurra il modello Hamman e Free-man(1977), che verra ricondotto alla sua forma base, cioe al modello “Competizione traSpecie ”. Si analizzeranno localmente tutti punti critici, definendone le caratteristichedi stabilita o instabilita. Quindi, al temine del nostro studio qualitativo, identificheremouna serie di condizioni che governano tutto il sistema ma, che allo stesso tempo, lascianoalcuni quesiti irrisolti. Per rispondere a queste domande costruiremo una Funzione diLyapunov per il problema in esame.

Indice

Introduzione 2

Elenco delle Figure. 7

1 Equazioni Differenziali non Lineari e Stabilita 81.1 Motivazioni . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 81.2 Piano delle Fasi: Sistemi Lineari . . . . . . . . . . . . . . . . . . . . . . . 91.3 Autovalori Reali Diversi Stesso Segno . . . . . . . . . . . . . . . . . . . . 101.4 Autovalori Reali di Segno Opposto . . . . . . . . . . . . . . . . . . . . . 121.5 Autovalori Reali Uguali . . . . . . . . . . . . . . . . . . . . . . . . . . . . 131.6 Autovalori Complessi . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 151.7 Autovalori Immaginari Puri . . . . . . . . . . . . . . . . . . . . . . . . . 171.8 Sommario Punti Critici . . . . . . . . . . . . . . . . . . . . . . . . . . . . 181.9 Sistemi Autonomi e Stabilita . . . . . . . . . . . . . . . . . . . . . . . . . 191.10 Metodo Indiretto di Lyapunov: Approssimazione di Sistemi non Lineari . 211.11 Perturbazioni di Sistemi Lineari . . . . . . . . . . . . . . . . . . . . . . . 251.12 Competizione tra Specie . . . . . . . . . . . . . . . . . . . . . . . . . . . 261.13 Modello Predatore-Preda . . . . . . . . . . . . . . . . . . . . . . . . . . . 35

2 Metodo Diretto di Lyapunov 402.1 Introduzione alla Teoria della Stabilita di Lyapunov . . . . . . . . . . . . 402.2 Teorema di Lyapunov . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 412.3 Costruzione della funzione di Lyapunov . . . . . . . . . . . . . . . . . . . 43

3 Descrizione e Significato del Modello Economico 463.1 Fenomeno Empirico . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 463.2 Modello del Fenomeno Empirico . . . . . . . . . . . . . . . . . . . . . . . 483.3 Analisi del Fenomeno Empirico: Ricerca dei Punti Critici . . . . . . . . . 503.4 Linearizzazione del Fenomeno Empirico . . . . . . . . . . . . . . . . . . . 523.5 Analisi Qualitativa: Fenomeno Empirico . . . . . . . . . . . . . . . . . . 553.6 Significato Economico dei Punti Critici . . . . . . . . . . . . . . . . . . . 58

5

3.7 Funzione di Lyapunov e Stabilita del Fenomeno Empirico . . . . . . . . . 59

Bibliografia 67

6

Elenco delle figure

1.1 Autovalori Reali Diversi Stesso Segno: Nodo . . . . . . . . . . . . . . . . 111.2 Autovalori Reali di Segno Opposto: Punto Sella . . . . . . . . . . . . . . 121.3 Due Autovettori Indipendenti: Punto Stella . . . . . . . . . . . . . . . . 141.4 Un Solo Autovettore Indipendente: Nodo Degenere . . . . . . . . . . . . 151.5 Autovalori Complessi Coniugati: Punto Spirale . . . . . . . . . . . . . . . 161.6 Autovalori Immaginari Puri: Centro . . . . . . . . . . . . . . . . . . . . . 171.7 Punti Critici . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 181.8 Quattro Casi Rilevanti per la Competizione tra Specie . . . . . . . . . . 321.9 4 Direzioni del Modello Predatore Preda . . . . . . . . . . . . . . . . . . 38

2.1 Interpretazione geometrica del metodo di Lyapunov . . . . . . . . . . . . 42

3.1 Schematizzazione Modello Mollona e Presutti(2006) . . . . . . . . . . . . 493.2 Risultati Matrice Hessiana in P12 . . . . . . . . . . . . . . . . . . . . . . 66

7

Capitolo 1

Equazioni Differenziali non Lineari eStabilita

1.1 Motivazioni

Le equazioni differenziali sono, uno dei piu importanti strumenti che l’analisi matematicamette a disposizione per lo studio di modelli applicati nei piu disparati settori scientifici:fisica, biologia, economia.Una comprensione di processi complessi e spesso ottenuta creando e combinando mod-elli base piu semplici. Quindi la conoscenza di questi modelli, delle equazioni e dellesoluzioni che li descrivono e il primo e indispensabile passo verso la soluzione di modelliad alta complessita che rappresentano problematiche reali. E’ dunque importante per larisoluzione di modelli prevalentemente complessi l’uso di una varieta di tecniche matem-atiche, analitiche e numeriche. Risultati quantitativi e grafici, prodotti dal computer,chepossono illustrare e chiarire conclusioni che potrebbero essere oscurate da espressionianalitiche particolarmente complesse.In altre parole, l’attuazione di procedure numeriche restera sempre un buon puntodi partenza per un’analisi preliminare; questa, sara la guida che ci portera verso ladefinizione di soluzioni qualitative e verso la scoperta di quei range di variabili che mer-itano grande importanza.Quindi per studiare tali modelli saranno richieste entrambe le tecniche analitiche e com-putative, sara compito di chi realizza lo studio scegliere quella maggiormente adatta alproblema in esame, ed utilizzarle entrambe per potenziare lo studio e la comprensionedel problema.

8

1.2 Piano delle Fasi: Sistemi Lineari

Il piano delle fasi associato a un sistema di equazioni differenziali, mostra il compor-tamento del sistema attraverso curve e linee piene, traiettorie del sistema al tempo t.La rappresentazione del piano delle fasi avviene attraverso i ritratti delle fasi, rapp-resentazioni geometriche delle traiettorie del sistema in cui ogni insieme di condizioniiniziali e rappresentato da una curva differente, o da un punto. Queste curve possonoessere pensate come spostamenti di una particella su diverse traiettorie la cui velocita edefinita dall’equazione differenziale. Tali ritratti sono di grandissima importanza, percherivelano informazioni non intuitive come traiettorie di convergenza(divergenza), bacinidi attrazione(repulsione) al variare del tempo. L’analisi e la tipicizzazione di questetraiettorie sara l’ingrediente base per l’utilizzo di una teoria qualitativa.

Data la difficolta computazionale di molte equazioni differenziali un approccio cal-colistico diretto non e sempre possibile e dunque diventa di grande importanza ottenereinformazioni qualitative senza attualmente risolvere le equazioni in esame.

Consideriamo l’equazione autonoma di primo grado

dy/dt = f(y), (1.1)

consideriamo inoltre con il sistema di equazioni differenziali piu semplice: omogeneodi secondo ordine a coefficienti costanti

dx/dt = Ax (1.2)

ove A e una matrice 2x2 e x un vettore 2x1. Sistemi di questo tipo vengono normal-mente risolti cercando soluzioni nella forma x = ξert. Sostituendo alla x dell’Eq (1.2)troveremo l’equazione.

(A− rI)ξ = 0. (1.3)

9

Un semplice calcolo mostra che r deve essere un autovalore mentre ξ e un autovet-tore della matrice A. Ricordiamo che gli autovalori della matrice A sono soluzionidell’equazione polinomiale

det(A− rI) = 0 (1.4)

mentre gli autovettori sono determinati risolvendo l’Eq.(1.3). Per l’Eq.(1.1) zero eun punto di grande importanza; questo punto corrisponde a una soluzione stazionariadetta soluzione di equilibrio o punto critico. Analogamente per il sistema (3.2),Ax = 0 corrisponde a una soluzione di equilibrio. Assumendo che A sia non-singolare eche detA 6= 0, l’unica possibile soluzione e x = 0.

Vogliamo ora studiare i diversi comportamenti del punto di equilibrio x = 0 a secondadelle tipologie degli autovalori reali distinti, reali coincidenti e complessi coniugati.

1.3 Autovalori Reali Diversi Stesso Segno

La soluzione generale per l’Eq.(1.2) dove r1 e r2 sono entrambi positivi o negativi e datada:

x = c1ξ(1)er1t + c2ξ

(2)er2t.

Supponiamo r1 < r2 < 0 e i corrispondenti autovettori ξ(1) e ξ(2) . Tutte le soluzionia prescindere dai valori assunti dalle costanti c1 e c2, si avvicinano al punto criticonell’origine per t → ∞. In questi casi e di grande utilita riscrivere l’equazione nelseguente modo.

x = er2t[c1ξ(1)e(r1−r2)t + c2ξ

(2)]

Scrivendo l’equazione generale in questa forma, si nota subito che mentre t → ∞,c1ξ

(1)er1t e trascurabile per t sufficientemente grande, in quanto r1 − r2 e minore di 0 equindi la soluzione seguira la traiettoria ξ(2) verso l’origine.

10

Tutte le altre soluzioni si avvicineranno all’origine tangenti a ξ(2), ad eccezione delletraiettorie che giacciono su ξ(1).

Figura 1.1: Autovalori Reali Diversi Stesso Segno: Nodo

Ora supponiamo analogamente che 0 < r2 < r1 e riscriviamo nuovamente la soluzionegenerale:

x = er1t[c1ξ(1) + c2ξ

(2)e(r2−r1)t],

si osserva che mentre t→∞ il termine c2ξ(2)e(r2−r1)t e trascurabile per t sufficiente-mente grande e che la soluzione segue la traiettoria ξ(1) allontanandosi dall’origine.Tutte le altre traiettorie si allontanano dall’origine e sono tangenti a ξ(1) ad eccezione diquelle che giacciono su ξ(2).

11

1.4 Autovalori Reali di Segno Opposto

La Soluzione generale per l’Eq.(1.2) e data da:

x = c1ξ(1)er1t + c2ξ

(2)er2t

ove stavolta supponiamo gli autovalori r1 > 0 e r2 < 0. Osserviamo il piano dellefasi.

Figura 1.2: Autovalori Reali di Segno Opposto: Punto Sella

Dal grafico vediamo che se la soluzione parte da un punto iniziale su ξ(1), c2 sarauguale a zero. Conseguentemente la soluzione rimarra su ξ(1) per qualunque t dato cher1 e maggiore di zero x → ∞ come t → ∞. Diversamente se la soluzione parte da unpunto iniziale su ξ(2), allora c1 sara uguale a zero visto che r2 e minore di zero, x → 0come t→∞.Quindi possiamo dire che le uniche traiettorie che si avvicineranno all’origine, sono quelleche giacciono precisamente su ξ(2), mentre tutte le altre seguiranno traiettorie simili aquelle mostrate in figura.

12

1.5 Autovalori Reali Uguali

Supponiamo r1 = r2 = r, consideriamo il caso in cui r1 e r2 siano entrambi nega-tivi(positivi traiettorie simili direzione inversa). Esistono due sottocasi: gli autovaloriripetuti hanno due corrispondenti autovettori linearmente indipendenti o uno solo.

• (a) Due autovettori indipendenti. La soluzione generale dell’Eq.(1.2) e nellaforma

x = c1ξ(1)er1t + c2ξ

(2)er2t

dove ξ(1) e ξ(2) sono autovettori indipendenti. Per tracciare le traiettorie e utilescrivere l’equazione in questa forma

x = ert[c1ξ(1) + c2ξ

(2)]

dalla quale possiamo osservare che il rapporto tra x2/x1 e indipendente da tma dipendente dagli autovalori e dalle costanti arbitrarie. L’esponenziale invecedefinisce la direzione delle traiettorie all’aumentare di t. Per tali motivi ognitraiettoria giace in una linea passante per l’origine come in figura.

13

Figura 1.3: Due Autovettori Indipendenti: Punto Stella

• (b) Un autovettore indipendente. La soluzione generale dell’Eq(1.2) in questocaso e

x = c1ξ(1)ert + c2(ξ te

rt + η ert)

dove ξ e l’autovettore e η e l’autovettore generalizzato associato a quello ripetuto.Per tracciare le traiettorie e ancora opportuno scrivere l’equazione in questa forma

x = ert[(c1ξ + c2η) + c2ξ t] = yert

dalla quale possiamo osservare, che l’esponenziale ha effetto solo sulla grandezzadi x, ma che insieme a r modifichera la direzione delle traiettorie verso l’interno ol’esterno. Notiamo anche che c2ξ t e dominante per t sufficientemente grande e chese c2 = 0 allora x = c1ξ e

rt.Procediamo disegnando le linee date da (c1ξ+c2η)+c2ξ t per osservarne la direzioneall’aumentare di t. Inoltre osserviamo che le suddette traiettorie quando t = 0passano per il punto c1ξ + c2η. Se t aumenta la traiettoria segue l’incremento dit ma insieme alla diminuzione della grandezza di x che quindi si avvicinera allozero a causa del decadimento dell’esponenziale ert quando r1 = r2 < 0. Se invece

14

r1 = r2 > 0, le traiettorie saranno similari ma la direzione questa volta sara versol’esterno.

Figura 1.4: Un Solo Autovettore Indipendente: Nodo Degenere

1.6 Autovalori Complessi

Supponiamo che gli autovalori siano del tipo λ± ıµ, dove λ e µ sono reali, λ 6= 0 e µ > 0.Questo tipo di sistemi e caratterizzati dalla seguente forma:

x′ =

(λ µ−µ λ

)x.

Ora introduciamo le traiettorie del sistema rappresentate in coordinate polari dateda r e θ, ricordiamo che r2 = x2

1 + x22 e tanθ = x2/x1.

r = ceλ t (1.5)

θ = −µ t+ θ0 (1.6)

15

Ogni sistema con autovalori complessi e λ 6= 0 e sempre una spirale intorno all’originela cui direzione verso l’interno o l’esterno varia a seconda del termine ceλ t .

Figura 1.5: Autovalori Complessi Coniugati: Punto Spirale

Dall’Eq.(1.5) se t→∞ con λ 6= 0 possiamo vedere che r → 0 se λ < 0 oppure r →∞se λ > 0. Le traiettorie quindi sono spirali che si avvicinano o allontanano dall’origine aseconda del segno di λ .Infine θ dall’Eq.(1.6) ci fornisce informazioni sul senso orario o antiorario del moto. Seµ > 0 mentre t → ∞, θ diminuira quindi il moto sara orario, viceversa se µ < 0 pert→∞ θ aumentera e il moto sara antiorario.

16

1.7 Autovalori Immaginari Puri

Se λ = 0 il nostro sistema assume la forma:

x′=

(0 µ−µ 0

)x

e ha autovalori del tipo ±ıµ. Utilizzando nuovamente le coordinate polari osserviamoche

r = c

θ = −µ t+ µ0

dove c e θ0 sono costanti. Da cui le traiettorie sono ellissi, centrate nell’origine conmovimento orario o antiorario a seconda che il termine µ sia > 0 oppure < 0. Questoparticolare punto critico viene chiamato centro.

Figura 1.6: Autovalori Immaginari Puri: Centro

17

1.8 Sommario Punti Critici

Figura 1.7: Punti Critici

18

1.9 Sistemi Autonomi e Stabilita

• Sistema Autonomo. Definiamo sistema autonomo un sistema composto da dueequazioni differenziali, nella seguente forma

{dx/dt = F (x, y)

dy/dt = G(x, y)

ove le funzioni F e G sono continue, con derivate parziali del primo ordine continuein un qualche dominio D del piano xy. Se (x0, y0) e un punto in D, allora esisteuna soluzione x = φ(t) , y = ψ(t) del sistema che soddisfa le condizioni iniziali

x(t0) = x0 , y(t0) = y0

definita in qualche intervallo aperto I di tempo contenente t0. Frequentemente ilproblema viene riscritto nella forma

dx/dt = f(x) , x(t0) =x 0

dove x= xi+yj, f(x) = F (x, y)i+G(x, y)j e x 0 = x0i+y0j e (i, j) la base canonicadi R2. La soluzione del sistema e espressa come x = Φ(t) dove Φ(t) = φ(t)i+ψ(t)j,che interpreteremo come una curva tracciata muovendo un punto sul piano xy.Osserviamo inoltre che F e G non dipendono direttamente dalla variabile t, masolamente attraverso x e y. Un esempio di sistema autonomo a due dimensioni edato dai sistemi lineari che assumono la forma:

d

dt

(xy

)=

(a bc d

)(xy

)

dove A e una matrice 2x2. In altre parole, se uno o piu elementi della matrice Asono in funzione della variabile indipendente t, il sistema sara non-autonomo.

19

• Stabilita e Instabilita. Sia

x′ = f(x)

I punti dove f(x) = 0 sono chiamati punti critici del sistema autonomo. In questipunti x′ = 0 e dunque questi punti corrispondono a soluzioni di equilibrio per ilsistema.Un punto critico x0, viene detto stabile se, data qualsiasi ε > 0, esiste una δ > 0tale che ogni soluzione x = Φ(t) del sistema, quando t = 0 soddisfa

‖Φ(0)− x0‖ < δ

ed inoltre soddisfa

‖Φ(t)− x0‖ < ε

per ogni t ≥ 0 in cui la soluzione esiste. Ogni soluzione che parte da un puntosufficientemente vicino (entro δ) a x0 restera vicino a x0(cioe entro la distanza ε).Se il punto critico non e stabile allora si dice instabile.

Un punto critico x0 si dice asintoticamente stabile se e stabile e se esiste δ0, con0 < δ0 < δ tale che se una soluzione soddisfa

‖Φ(0)− x0‖ < δ0

allora

limt→∞

Φ(t) = x0

Le traiettorie vicine a x0 non restano solamente in prossimita, ma si avvicinano ax0 per t→∞.

20

1.10 Metodo Indiretto di Lyapunov: Approssimazione

di Sistemi non Lineari

Si consideri il sistema non lineare autonomo in due dimensioni

x′ = f(x). (1.7)

In primis vogliamo descrivere le traiettorie del sistema nelle vicinanze del punto criticox0. Per farlo approssimeremo il sistema non lineare con uno lineare piu semplice daanalizzare e descrivere.

E sempre conveniente scegliere un punto x0 6= 0 vicino all’origine. Sara poi semprepossibile sostituire al punto x = x0 + u nell’Eq.(1.7), soddisfando il sistema nuovamentecon un punto nell’origine.

Si consideri un sistema non lineare

x′ = Ax+ g(x) (1.8)

tale che x = 0 sia un punto critico isolato per il sistema (1.7). Questo significa cheesiste un qualche intorno del punto, all’interno del quale non esistono altri punti critici.Assumiamo anche che det(A) 6= 0 cosı che x = 0 sia punto critico anche per il sistemalineare x′ = Ax.Per il sistema non lineare, avvicinarsi al sistema lineare, significa che g(x) e molto piccoloo piu precisamente g(x) ha derivate parziali continue del primo ordine e soddisfa lacondizione sui limiti.

‖g(x)‖/‖x‖ → 0 x→ 0

Si osserva che ‖g(x)‖ e piccolo rispetto a ‖x‖ vicino all’origine. Questo tipo di sistemae chiamato quasi lineare in prossimita di x = 0.

21

ESEMPIO

(xy

)′=

(1 00 0.5

)+

(xy

)+

(−x2 − xy

−0.75xy − 0.25y2

)

Un punto critico e (0, 0) e il det(A) 6= 0. Inoltre non e difficile determinare che glialtri punti critici sono (0, 2), (1, 0) e (0.5, 0.5), conseguentemente l’origine e un puntocritico isolato. Si verifichi che g(x) rispetti la condizione sui limiti

g1(x, y)

r=−x2 − xy

r=−r2 cos2 θ − r2 cos θ sin θ

r= −r(cos2 θ + cos θ sin θ)→ 0

come r → 0. Allo stesso modo anche g2(x,y)r→ 0 come r → 0 quindi il sistema e

quasi-lineare.

22

Riscriviamo il sistema di Eq.(1.7) come:

x′ = F (x, y) , y′ = G(x, y)

Se assumiamo che le funzioni F e G abbiano derivate parziali continue fino al secondoordine, possiamo dare la seguente approssimazione lineare del sistema, per un intorno diun punto critico (x0, y0). Consideriamo lo sviluppo in serie di Taylor nel punto criticoper le funzioni F e G.

F (x, y) = F (x0, y0) + Fx(x0, y0)(x− x0) + Fy(x0, y0)(y − y0) + η1(x, y)

G(x, y) = G(x0, y0) +Gx(x0, y0)(x− x0) +Gy(x0, y0)(y − y0) + η2(x, y)

dove η1(x, y)/[(x − x0)2 + (y − y0)

2]1/2 → 0 come (x, y) → (x0, y0) e analogamenteper η2. Notiamo inoltre che F (x0, y0) = G(x0, y0) = 0 e che dx/dt = d(x − x0)/dt edy/dt = d(y − y0)/dt. Dunque in un intorno del punto critico (x0, y0) il nostro sistemae:

d

dt

(x− x0

y − y0

)=

(Fx(x0, y0) Fy(x0, y0)Gx(x0, y0) Gy(x0, y0)

)(x− x0

t− y0

)+

(η1(x, y)η2(x, y)

)

Il sistema lineare che approssima il sistema non-lineare in prossimita di (x0, y0), edefinito dalla parte lineare dell’Eq.(1.11):

d

dt

(uv

)=

(Fx(x0, y0) Fy(x0, y0)Gx(x0, y0) Gy(x0, y0

)(uv

)(1.9)

dove u = x − x0 e v = y − y0. Vediamo un esempio pratico di come sia possibileapprossimare un sistema non lineare in un intorno di un suo punto critico.

23

ESEMPIO

Consideriamo il sistema non lineare:

{dx/dt = F (x, y) = y

dy/dt = G(x, y) = −ω2 sin x− γ y

I punti critici sono (0, 0) e (π, 0). Poiche:

Fx = 0 Fy = 1 Gx = −ω2 cos x Gy = −γ

nell’origine il corrispondente sistema lineare che approssima il sistema non-lineare e:

d

dt

(xy

)=

(0 1−ω2 −γ

)(xy

)

Analogamente per il punto (π, 0) il sistema lineare corrispondente e dato da:

d

dt

(uv

)=

(0 1ω2 −γ

)(uv

)(1.10)

dove x = π + u , y = v. E’ ragionevole pensare che almeno vicino all’origine cioche trasciniamo sia piccolo comparato ad Ax e quindi il nostro sistema sia una buonaapprossimazione.Sostanzialmente possiamo dire che per piccole variazioni della x, il sistema non lineareha a sua volta piccole perturbazioni che di fatto non hanno effetto sulle caratteristichedel punto critico e sulla stabilita del sistema, ad eccezione dei due casi sensibili. Ossiar1 e r2 sono puri immaginari oppure r1 e r2 sono reali e coincidenti.E’ quindi ragionevole aspettarsi che solo in questi due casi piccole perturbazioni alterinola stabilita e che in tutti gli altri sia possibile studiare il sistema attraverso un semplicesistema lineare.

24

1.11 Perturbazioni di Sistemi Lineari

Dato il sistema linearex′ = Ax (1.11)

e ricordando che se il det(A) 6= 0, x = 0 sara il solo punto critico per il sistema.Il punto critico si dira asintoticamente stabile se r1 e r2 sono valori reali negativi ohanno una parte reale negativa; stabile ma non asintoticamente stabile se se r1 e r2 sonopuri immaginari; instabile se r1 e r2 sono reali e l’uno o l’altro sono positivi o hannouna parte reale positiva.

E’ ormai chiaro che gli autovalori r1 e r2 della matrice A determinano il tipo di puntocritico in x = 0 le sue traiettorie e conseguentemente le su caratteristiche di stabilita oinstabilita.

Questi sistemi spesso vengono utilizzati in campi applicati e quindi con misure fisichea meno di una certa incertezza. Diventa quindi importante comprendere come piccolevariazioni perturbino le caratteristiche di stabilita o instabilita del punto critico.In generale piccole modifiche dei coefficienti si traducono in altrettanto piccole pertur-bazioni degli eigenvalues, strade dell’equazione polinomiale det(A− rI) = 0.

Esistono tuttavia eccezioni, situazioni sensibili come quando r1 = iµ e r2 = −iµ. Ilpunto critico e il centro, mentre le traiettorie sono curve chiuse intorno al punto critico.Dopo una piccola variazione r′1 = λ

′1 + iµ′ e r′2 = λ

′1 − iµ′, se λ

′1 e debole allora µ ∼= µ′,

altrimenti se λ′ 6= 0 come spesso accade e se λ′ < 0 il sistema puntera il centro restandostabile, altrimenti se λ′ > 0 instabile.Un altro caso si presenta quando r1 = r2; in questa situazioni piccole perturbazioni neicoefficienti generalmente causano una biforcazione delle soluzioni.Se le nuove soluzioni sono reali allora il punto critico del sistema perturbato restera unnodo ma se le soluzioni diventano complessi coniugati il punto critico diventera un nodospirale.

25

1.12 Competizione tra Specie

In questo capitolo vogliamo esplorare ed analizzare le interazioni legate all’analisi dialcuni problemi relativi a popolazioni dinamiche. Questi problemi sorgono quando dueo piu popolazioni interagiscono tra loro.

Siano x e y le popolazioni corrispondenti a due specie in competizione, per un certonumero limitato di scorte di cibo al tempo t. Supponiamo che ogni specie in assenza dirivali sia governata da una equazione.

{dx/dt = x(ε1 − σ1x)

dy/dt = y(ε2 − σ2y)

Con ε1 e ε2 tassi di crescita della popolazione e ε1/σ1 , ε2/σ2 livelli di saturazione.Per mettere in competizione le specie bastera modificare ε1 − σ1x con ε1 − σ1x − α1ydove α1 e una stima del grado di interferenza della specie contrapposta. Otteniamo cosıil sistema:

{dx/dt = x(ε1 − σ1x− α1y)

dy/dt = y(ε2 − σ2y − α2x)

Poiche questo sistema modella una situazione concreta, ci interesseremo soltanto allesoluzioni x e y non negativi. Questo particolare modello sara alla base della nostrasperimentazione. Inizieremo a trattarlo partendo da un caso specifico, spostandoci poiverso la soluzione in forma parametrica, di grande utilita per la risoluzione di modelli dimaggiore complessita che tratteremo nell’ultimo capitolo.

26

ESEMPIO DI COMPETIZIONE TRA SPECIE

Consideriamo il sistema:

{dx/dt = x(1− x− y)

dy/dt = y(0.5− 0.25y − 0.75x)

I punti critici sono determinati ponendo f(x) = 0 da cui

{x(1− x− y) = 0

y(0.5− 0.25y − 0.75x)

Dopo un breve calcolo si ottiene che punti sono (0, 0), (1, 0), (0, 2) e (0.5, 0.5), i qualicorrispondono ai punti di equilibrio del nostro sistema. Si noti che i primi tre punticorrispondono all’estinzione di una o entrambe le specie e che il punto (0.5, 0.5) e l’unicoin cui coesistono le due specie. La sopravvivenza delle specie, e quindi determinata dallostato iniziale del sistema. Per analizzare la stabilita studiamo l’approssimazione linearedel sistema nelle vicinanze di ogni punto critico.

Punto critico x = 0, y = 0.

d

dt

(xy

)=

(1 00 0.5

)(xy

)

Gli autovalori sono determinati risolvendo:

(1− r 0

0 1/2− r

)= 0 ,

27

cioe

r2 − 3/2r + 1/2 = 0

r1,2 =3/2±

√1/4

2

r1 = 1 r2 = 0.5 .

I corrispondenti autovettori sono dati da

(0 00 −1/2

)(ξ1ξ2

)=

(00

),

(1/2 00 0

)(ξ1ξ2

)=

(00

)

Da cui risulta subito che :

ξ(1) =

(10

)e ξ(2) =

(01

).

Quindi la soluzione generale e

(xy

)= c1

(10

)et + c2

(01

)e1/2t

L’origine e una soluzione instabile per il sistema lineare e conseguentemente ancheper il sistema non lineare. Tutte le traiettorie lasciano l’origine tangenti all’asse y, adeccezione della traiettoria che giace sull’asse x.

28

Punto critico x = 1, y = 0.

Sostituendo all’equazione x = 1 +u e y = v il corrispondente sistema lineare traslatosara:

d

dt

(uv

)=

(−1 −10 −0.25

)(uv

)

I suoi autovalori e autovettori sono

r1 = −1, ξ(1) =

(10

), r2 = −0.25, ξ(2) =

(4−3

),

e la sua soluzione generale e data da:

(uv

)= c1

(10

)e−t + c2

(4−3

)e−0.25t .

Il punto (1, 0) e asintoticamente stabile per il sistema lineare e conseguentementeanche per il sistema non lineare. Se i valori iniziali di x e y sono sufficientemente vicinial punto, il processo di iterazione portera all’estinzione della specie y. Osserviamo cheesistono un paio di traiettorie che approcciano il punto sull’asse x mentre tutte le altresi avvicineranno a (1, 0) tangenti a ξ(2).

Punto critico x=0 , y=2.

Procedendo in modo analogo a quanto fatto precedentemente otteniamo:

d

dt

(uv

)=

(−1 0−1.5 −0.5

)(uv

)

I suoi autovalori e autovettori sono dati da

r1 = −1, ξ(1) =

(13

), r2 = −0.5, ξ(2) =

(01

),

e la sua soluzione generale e

29

(uv

)= c1

(13

)e−t + c2

(01

)e−0.5t ,

Il punto critico (0, 2) e asintoticamente stabile come il punto precedente. Tutte letraiettorie si avvicinano al punto critico lungo l’asse y ad eccezione della traiettoria chesi avvicina lungo l’autovettore ξ(1).

Punto critico x=0.5 , y=0.5.

Ragionando come in precedenza abbiamo:

d

dt

(uv

)=

(−0.5 −0.5−0.375 −0.125

)(uv

).

Gli autovalori e autovettori sono:

r1 = −5+√

5716

∼= 0.1594, ξ(1) =

(1

(−3−√

57)/8

),

r2 = −5−√

5716

∼= −0.7844, ξ(2) =

(1

(−3 +√

57)/8

),

e la soluzione generale e data da:

(uv

)= c1

(1

(−3−√

57)/8

)e0.159t + c2

(1

(−3 +√

57)/8

)e−0.7844t

Poiche gli autovalori hanno segno opposto, il punto critico e un punto sella, quindiinstabile. Un paio di traiettorie si avvicinano al punto critico per t → ∞, mentre lerestanti si allontanano da esso. Si osserva inoltre che mentre le traiettorie entranti siavvicinano al punto critico, sono tangenti all’autovettore ξ(1).

30

ESEMPIO COMPETIZIONE TRA SPECIE SISTEMA GENERICO

L’esempio precedente mostra chiaramente come in alcuni casi la competizione traspecie approdi verso uno stato di coesistenza, mentre in altri la competizione porti al-l’estinzione di una delle specie.Per comprendere pienamente il meccanismo che porta a uno scenario piuttosto che all’al-tro, e utile analizzare nella forma piu generale possibile le funzioni F e G che definisconol’interazione.Quattro sono i casi rilevanti da tenere in considerazione, dipendenti dalla direzione dellerette.

x(ε1 − σ1x− α1y) = 0

y(ε2 − σ2y − α2x) = 0

Risolvendo l’equazione algebrica osserviamo che l’uguaglianza e verificata da x = 0 ex = (ε1 − α1y)/σ1. Risolviamo ogni caso singolarmente.

Caso (x = 0)

y(ε2 − σ2y) = 0

L’uguaglianza e verificata per y = 0 e y = (ε2)/σ2, i primi due punti critici sono

P1 = (0, 0) , P2 = (0, ε2σ2

)

31

Caso (x = ε1−α1yσ1

)

yε2 − σ2y2 − α2

σ1

y (ε1 − α1) = 0

y( y (−σ1σ2 + α1α2) + (σ1ε2 − α2ε1) ) = 0

l’uguaglianza e verificata per y = 0 e y = (σ1ε2 − α2ε1)/(σ1σ2− α1α2), che sostituitialla prima equazione ci danno x = ε1/σ1 e x = (σ2ε1 − α1ε2)/(σ1σ2 − α1α2), gli ultimidue punti critici sono.

P3 = ( ε1σ1, 0) , P4 = ( σ2ε1−α1ε2

σ1σ2−α1α2, σ1ε2−α2ε1σ1σ2−α1α2

)

Figura 1.8: Quattro Casi Rilevanti per la Competizione tra Specie

32

Supponiamo che (X, Y ) sia uno qualunque dei quattro punti critici precedentementedeterminati con X e Y entrambi positivi. Per studiare il sistema nelle vicinanze di questipunti dobbiamo osservare il sistema lineare ottenuto da quello non lineare.

Fx(X, Y ) = ε1 − 2σ1x− αy , Fy(X, Y ) = −α1X

Gx(X, Y ) = −α2Y , Gy(X, Y ) = ε2 − 2σ2Y − α2X

Il sistema lineare che approssima del sistema non lineare e data da:

d

dt

(uv

)=

(ε1 − 2σ1x− αy −α1X−α2Y ε2 − 2σ2Y − α2X

)(uv

).

L’equazione. determinata stabilisce le condizioni sotto le quali e possibile la coesisten-za delle specie x e y. Dei quattro casi, solo il (c) e (d) possono portare a questo stato,cioe solo nel punto critico P4 sara possibile la coesistenza nel lungo periodo.

X =σ2ε1 − α1ε2σ1σ2 − α1α2

, Y =σ1ε2 − α2ε1σ1σ2− α1α2

Nel punto critico ε1 − σ1X − α1Y = 0 e ε2 − σ2Y − α2X = 0 e quindi il sistemalineare generale potra essere ridotto alla forma.

d

dt

(uv

)=

(−σ1X −α1X−α2Y −σ2Y

)(uv

)

33

Risolvendo il polinomio caratteristico otteniamo gli eigenvalues del sistema.

det(

(−σ1X − r −α1X−α2Y −σ2Y − r

)) = 0

r2 + r(σ1X + σ2Y ) + (XY )(σ1σ2 − α1α2) = 0

r1,2 =−b±

√b2 − 4ac

2a=−(σ1X + σ2Y )±

√(σ1X + σ2Y )2 − 4(XY )(σ1σ2 − α1α2)

2

Se σ1σ2 − α1α2 < 0, il radicando e positivo e maggiore di (σ1X + σ2Y )2. Da cui gliautovalori saranno reali e di segno opposto, conseguentemente il punto sara un puntosella, la coesistenza sara impossibile. Se σ1σ2 − α1α2 > 0, il radicando sara minore di(σ1X + σ2Y )2. In questo caso gli autovalori saranno reali, negativi e diversi tra lorooppure complessi con una parte reale negativa. Un attenta analisi inoltre del radicandomostra che valori complessi non sono ottenibili quindi il punto critico e asintoticamentestabile e quindi la coesistenza e possibile.Ricordando la condizione che vuole X, Y positivi, quando σ1σ2 > α1α2, la competizionee debole e le specie possono coesistere; quando σ1σ2 < α1α2, l’interazione e forte e quindile specie non potranno coesistere, una dovra morire.

Supponendo stabile il centro del sistema, si osserva che per gli altri punti dovrannovalere le seguenti condizioni.

ε1σ2 − α1ε2 > 0 e σ1ε2 − α2ε1 > 0

Tale per cui i punti critici qualitativamente sono.

• il punto (0, 0) e instabile

• il punto ( ε1σ1, 0) e instabile

• il punto (0, ε2σ2

) e instabile

34

1.13 Modello Predatore-Preda

Questo modello rappresenta una evoluzione del modello precedente. Il modello predatorepreda sotto analisi riprendera un situazione di interazione tra due specie; una (il preda-tore), preda l’altra (la preda). Le due specie sono sostentate da due diverse tipologiedi cibo. Chiaramente un modello cosı semplice non puo ricreare un fenomeno di talecomplessita, tuttavia puo essere un primo passo verso una maggiore comprensione delfenomeno. Siano x e y popolazioni di prede e predatori al tempo t.

Assumiamo che :

• In assenza di predatori le prede aumentino in percentuale proporzionalmente allivello corrente della popolazione: dx/dt = ax , a > 0 quando y = 0

• In assenza di prede i predatori muoiano: dy/dt = −cy , c > 0 quando x = 0

• Il numero degli incontri tra le due specie e proporzionale al prodotto delle loropopolazioni. Ogni incontro tende ad aumentare i predatori e diminuire le prede.Quindi l’incremento dei predatori e di γ xy mentre le prede diminuiscono di −αxycon α e γ costanti positive.

Il sistema e dato da:

{dx/dt = x(a− α y)

dy/dt = y(−c+ γ x)

anche detta Eq. di Lotka-Volterra.

35

Per calcolare i punti critici poniamo f(x) = 0

{x(a− α y) = 0

y(−c+ γ x) = 0

I punti critici sono (0, 0) e (c/γ, a/α). Si valuti la soluzione del corrispondente sistemalineare in ogni punto critico.

Punto critico (0,0)

d

dt

(xy

)=

(a 00 −c

)(xy

)

La matrice A e una matrice diagonale e quindi gli autovalori saranno esattamente icoefficienti della diagonale. Gli autovalori e i corrispondenti autovettori sono

r1 = a , ξ(1) =

(10

); r2 = −c , ξ(2) =

(01

)

e la soluzione generale e

(xy

)= c1

(10

)eat + c2

(01

)e−ct .

Si osserva che a > 0 mentre −c < 0 con a, c costanti positive. Ne segue che il puntocritico (0, 0) e un punto sella per il sistema lineare e quindi anche per il corrispondentesistema non lineare, in cui solo due traiettorie si avvicinano all’origine sull’asse y, mentretutte le altre si allontanano.

36

Punto critico ( cγ, aα

)

Sostituendo x = (c/γ) + u, y = (a/α) + v troviamo il corrispondente sistema lineare:

d

dt

(uv

)=

(0 −α c

γ

γ aα

0

)(uv

)

e che gli autovalori e i corrispondenti autovettori sono immaginari puri, della formar = ± i

√ac. Il punto critico e quindi stabile e un centro per il sistema. Per identificare le

traiettorie e possibile dividere la seconda equazione del sistema per la prima. Ottenendo

dv

du=dv/dt

du/dt=

γ aαu

−α cγv

α2cv dv + γ2au du = 0

infine integrando:

α2c v2 + γ2a u2 = k

con k costante di integrazione. Le traiettorie che circondano il punto critico ( cγ, aα

),sono ellissi ed e possibile dimostrare che per k costante il punto critico e anche un centroper il sistema non lineare.

37

Figura 1.9: 4 Direzioni del Modello Predatore Preda

La variazione ciclica puo essere analizzata nel dettaglio quando deviazioni dal puntocritico sono piccole e quindi

u =c

γK cos(

√act+ φ)

v =a

α

√c

aK sin(

√act+ φ)

e

x =c

γ+c

γK cos(

√act+ φ)

y =a

α+a

α

√c

aK sin(

√act+ φ)

38

Queste equazioni sono buone approssimazioni per le traiettorie ellittiche nell’intornodel centro del sistema. Osservandole vediamo che:

• Il numero dei predatori e delle prede varia con un periodo di 2π/√ac indipenden-

temente dalle condizioni iniziali.

• I predatori e le prede sono fuori fase per un quarto di ciclo. Le prede conducono ei predatori le seguono.

• L’ampiezza dell’oscillazione e di Kc/γ per le prede e a√cK/α

√a per i predatori.

Inoltre dipende dalle condizioni e dai parametri del problema.

• La popolazione media dopo un ciclo e la stessa presente nel punto di equilibrio( cγ, aα

).

Il modello predatore preda di Lotka-Volterra dimostra come sia possibile antici-pare variazioni cicliche. Attraverso il suo studio e possibile giungere a conclusioni nonintuitive, applicabili nei piu svariati campi.

39

Capitolo 2

Metodo Diretto di Lyapunov

2.1 Introduzione alla Teoria della Stabilita di Lya-

punov

Il problema della stabilita di un sistema dinamico e una questione di fondamentale im-portanza, sia teorica che pratica. Come abbiamo piu volte ripetuto un sistema e stabilese la sua dinamica, a partire da certe condizioni iniziali nell’intorno di un punto criticoper il sistema, raggiunge tale punto e in questo caso possiamo dire che non si evolve ulte-riormente. Un importante approccio alla stabilita di sistemi non lineari, fu introdotto daAlexander Mikhailovich Lyapunov alla fine del secolo diciannovesimo. L’approcciodi Lyapunov al problema comprende due metodi: il metodo indiretto che abbiamo giatrattato (linearizzazione) e il metodo diretto.

Il metodo indiretto determina le caratteristiche di stabilita locale di un sistema nonlineare nell’intorno di un punto di equilibrio di un sistema non lineare, attraverso lo studiodella sua approssimazione lineare vicino a tale punto. Sfortunatamente questo metodo euna teoria locale e quindi non e in grado di determinare il bacino di attrazione del puntocritico, cioe un dominio contenente il punto stesso e tale che tutte le traiettorie del sistemache si originano all’interno di tale dominio e tendano ad avvicinarsi ed eventualmenteraggiungere il punto critico per t→∞.

Il metodo diretto e nato proprio per affrontare il problema della determinazione delbacino di attrazione. Questo metodo affronta il problema della stabilita globale deisistemi non lineari, generalizzando due principi della fisica. Lyapunov infatti formalizzol’idea secondo cui se l’energia totale di un sistema viene dissipata il sistema deve esserestabile. Il vantaggio del metodo diretto e dato dal fatto che le caratteristiche di stabilita

40

di un sistema non lineare di complessita anche molto elevata possono essere definitestudiando una funzione ausiliaria scalare del tipo energia, senza quindi dover risolvereanaliticamente le equazioni differenziali che descrivono la dinamica del sistema.

Tale metodo e molto potente e ci garantisce un gran numero di informazioni comeper esempio l’estensione della base di attrazione anche se non e sempre facile da attuare.

2.2 Teorema di Lyapunov

Nel metodo diretto, come vedremo, la funzione di Lyapunov riveste un importanzafondamentale. Vediamo le sue caratteristiche e i tre teoremi fondamentali di Lyapunov.

Funzione Definita Positiva/Negativa

Sia V : D → R una funzione definita in un dominio D ⊂ Rn, ove D contiene l’origine.Diremo che V e definita positiva in D se V (0, 0) = 0 e V (x, y) > 0 per tutti gli altripunti in D.Similarmente diremo che V e definita negativa se V (0, 0) = 0 e V (x, y) < 0 perogni altro punto del dominio. Infine se la disuguaglianza e non stretta (≤≥) allorachiameremo la funzione semi definita positiva o semi definita negativa.

ESEMPIO FUNZIONI DEFINITE POSITIVE/NEGATIVE

Consideriamo la funzione energia V : R2 → R

V (x, y) = sin(x2 + y2)

Vediamo subito che essa e definita positiva nel dominio 0 < x2 + y2 < π/2 conV (0, 0) = 0, infatti in tale dominio V (x, y) > 0. Vediamo la funzione:

V (x, y) = (x+ y)2 .

Questa funzione e solamente semi-definita positiva infatti V (0, 0) = 0 e y = −x seV (x, y) = 0.

41

Primo Teorema di Lyapunov

Supponiamo che il sistema autonomo abbia un punto critico isolato all’origine. Se esisteuna funzione V continua detta funzione energia con derivate parziali del primo ordinecontinue positivamente definite e per le quali V (x, y) definita come Vx(x, y)dx

dt+Vy(x, y)dy

dt

e negativamente definita in un qualche dominio D nel piano-xy contenente (0,0), alloral’origine e asintoticamente stabile. Se V e negativa ma solamente semi-definita alloral’origine e un punto critico stabile.

Figura 2.1: Interpretazione geometrica del metodo di Lyapunov

Secondo Teorema di Lyapunov

Supponiamo che l’origine sia un punto critico isolato per il sistema autonomo conside-riamo V funzione continua con derivate parziali continue V : D → r. Si presume cheV (0, 0) = 0 e che nelle vicinanze dell’origine ci sia almeno un punto per il quale V epositiva (negativa). Se esiste un dominio D contenente l’origine tale che la funzione V epositiva definita(negativa definita) in D, allora l’origine e un punto critico instabile.

Terzo Teorema di Lyapunov

Supponiamo che l’origine sia un punto critico isolato per il sistema autonomo e conside-riamo V funzione continua con derivate parziali continue. Se esiste un dominio limitato

42

Dk contenente l’origine dove V (x, y) < K, V e positiva definita e V e negativa defini-ta, allora ogni soluzione del sistema omogeneo che parte da un punto in Dk si avvicinaall’origine per t che tende all’infinito.

2.3 Costruzione della funzione di Lyapunov

La difficolta maggiore del il metodo diretto di Lyapunov e che non esiste un metodo peridentificare la funzione scalare V che ci permette di analizzare il sistema attraverso i treteoremi di Lyapunov.Nel caso in cui il sistema autonomo rappresenti un problema della fisica e naturaleconsiderare come funziono V , la funzione dell’energia totale.Sfortunatamente in tutti gli altri casi non possiamo avvalerci di alcun metodo generaleper costruire la funzione; per alcune classi speciali di equazioni ci si e resi conto dell’utilitadel seguente teorema per creare funzioni definite negative o positive.

La funzione

V (x, y) = ax2 + bxy + cy2

• e positiva definita se e solo se a > 0 e 4ac− b2 > 0

• e negativa definita se e solo se a < 0 e 4ac− b2 < 0

43

Esempio del Metodo Diretto: Competizione tra Specie

dx/dt = x(1− x− y)

dy/dt = y(0.75− y − 0.5x)

Riprendendo i concetti precedentemente trattati, questo sistema modella un certotipo di competizione tra specie che nel punto critico (0.5,0.5) il sistema e asintotica-mente stabile. Utilizzando una funzione di Lyapunov vogliamo determinare un bacinodi attrazione per tale punto.

Sostituendo x = 0.5 + u , y = 0.5 + v alle equazione del sistema omogeneo.

du/dt = −0.5u− 0.5v − u2 − uv

dv/dt = −0.24u− 0.5v − 0.5uv − v2

Si consideri come funzione di Lyapunov V (u, v) = u2 + v2 che chiaramente e definitapositiva e in V (0, 0) = 0. Ora dobbiamo solamente determinare una regione contenentel’origine nel piano-uv, nella quale le derivate V rispetto al sistema siano definite negative.

V (u, v) =

2u(du/dt) + 2v(dv/dt)

− [(u2 + 1.5uv + v2) + (2u3 + 2u2v + uv2 + 2v3)]

44

Vogliamo mostrare che l’espressione tra parentesi quadre e definita positiva almenoper u, v sufficientemente piccoli.Se scriviamo u2 + 1.5uv + v2 = 0.25(u2 + v2) + 0.75(u+ v)2, vediamo che questi terminisono definiti positivi.Ora cerchiamo di mostrare che nelle vicinanze di u = 0, v = 0 i termini cubici sono piupiccoli in forza dei termini quadrati |2u3 +2u2v+uv2 +2v3| < 0.25(u2 +v2)+0.75(u+v)2.

|2u3 + 2u2v + uv2 + 2v3| = r3|2 cos3 θ + 2 cos2 θ sin θ + cos θ sin2 θ + 2 sin3 θ|

|2u3 + 2u2v + uv2 + 2v3| = 7r3

Dato che | sin θ|, | cos θ| sono ≤ 1 diremo che l’equazione e definita positiva se everificata la disuguaglianza.

7r3 < 0.25(u2 + v2) = 0.25r2

r < 1/28

L’origine e quindi asintoticamente stabile e di conseguenza anche il punto critico(0.5, 0.5) quando r < 1/28, cioe abbiamod eterminato un bacino di attrazione circolaredi raggio r = 1/28 centrato in 0.5, 0.5).

45

Capitolo 3

Descrizione e Significato del ModelloEconomico

3.1 Fenomeno Empirico

Vogliamo prendere in considerazione la realta di un distretto industriale, e il suo sviluppoe dare una spiegazione alla dinamica evolutiva delle popolazioni delle organizzazionilocalizzate, ossia quelle aggregazioni che subiscono il fenomeno dell’internazionalizzazionepassiva.In altre parole prenderemo in considerazione il fenomeno della competizione tra impresee il loro cambiamento a seguito dell’ingresso di nuovi competitori stranieri.

Nel lavoro (Mollona e Presutti, 2006) gli autori esaminano un distretto industrialetessile italiano. Comincino la loro osservazione da una situazione di equilibrio interno,che nel corso degli anni 90 si e via via destabilizzata, a causa della localizzazione velocedi piccoli fornitori specializzati stranieri, procedendo verso una concorrenza forte frafornitori locali e fornitori stranieri in un’area limitata.

Negli anni 90 infatti molti distretti italiani hanno dislocato la loro produzione all’es-tero internazionalizzandosi ossia aprendosi a un processo di frammentazione territorialedella catena di valore.Questo processo ha comportato a sua volta una forte transizione del controllo internodel distretto verso di un piu aperto controllo globalmente integrato. Improvvisamentesono diventate rilevanti nuove strategie per migliorare il controllo sui mercati interni estranieri in un ottica di apertura internazionale.

46

In questo processo di globalizzazione, si osserva la tendenza all’aumento dei competi-tori stranieri. Questo nuovo sviluppo e stato chiamato internazionalizzazione passiva deidistretti italiani e come effetto ha prodotto un conflitto potenziale fra la rete locale ormaiconsolidata e i nuovi entrati stranieri.

In generale le domande alle quali vorremmo trovare una risposta sono le seguenti:

• Come possiamo valutare l’entrata di diversi attori stranieri nei distretti italiani?

• Per quale motivo gli attori stranieri erano in grado di sostituire velocemente glioperatori locali tradizionali?

• Qual’e la dinamica e l’evoluzione della concorrenza fra le popolazioni delle ditteesistenti e le popolazioni dei nuovi arrivati?

• Quali modifiche ha portato nel modello di controllo tradizionale del distretto l’en-trata di aziende straniere?

La dinamica di questo processo e chiaramente ambigua e di difficile giudizio datala carenza di dati empirici che mettano in evidenza l’equilibrio dei distretti sotto l’ef-fetto dell’internazionalizzazione passiva. Sin dagli anni 60 l’industria tessile italiana hasviluppato caratteristiche specifiche di qualita e manualita che hanno determinato unparticolare modello industriale basato sulla presenza di piccole e medie imprese special-izzate.E’ quindi interessante considerare l’ingresso di nuove popolazioni nel distretto industrialee cercare di spiegarne le interazioni.

47

3.2 Modello del Fenomeno Empirico

Riprendendo il fenomeno empirico precedentemente descritto della “Competizione traSpecie ”, cercheremo di definire un modello che descriva la concorrenza fra le popolazionidei fornitori e dei produttori delle merci all’interno di una serie di ingranaggi industriali,osservando il fenomeno da una prospettiva di ecosistema biologico.

Nell’ambiente in esame all’interno di un distretto industriale si osserva che i produttoridi beni finiti generano risorse per la popolazione dei fornitori.Consideriamo i fornitori tradizionali originari del distretto industriale ed i nuovi fornitoriesterni entranti come popolazioni eterogenee, dato lo sviluppo relazionale e le strategiecompetitive differenti che comportano reazioni diverse alle dinamiche ambientali esogene.

D’altra parte, anche i produttori delle merci finite che appartengono al distretto equelli che arrivano dall’esterno, sono due popolazioni differenti dalle organizzazioni condifferenti costi e prezzi di vendita, ma che a loro volta interferiscono nello stesso distrettoindustriale delle risorse creando competizione.

Il modello “Competizione tra Specie ”, (Hamman e Freeman(1977)), e definito daequazioni nella forma:

dxidt

= gixi(u− cji − xi

u)

Nell’equazione, gli indici i e j rappresentano le popolazioni in competizione. gi e iltasso di crescita della popolazione i mentre u e la risorsa limitata disponibile nel distrettonecessaria alla sopravvivenza degli attori.Infine cji e il tasso di competitivita della popolazione j rispetto alla popolazione i ossiala probabilita che un membro della popolazione j batta un membro della popolazione inell’acquisizione delle risorse.

48

Utilizzando come base il modello Hamman e Freeman(1997), Mollona e Presutti(2006)ne propongono una ritrattazione ampliata in cui la competizione si svolge su due livelli.

Piu precisamente, la quantita u e il numero di imprese clienti delle popolazioni all’in-terno dello stesso distretto. Conseguentemente, denominando con x1, x2 le popolazionidei fornitori appartenenti allo stesso distretto, con y1, y2 popolazioni di produttori dibeni finiti dello stesso distretto ed u numero di utenti finali potenziali acquirenti delprodotto, avremo un modello competitivo schematizzato in figura.

Figura 3.1: Schematizzazione Modello Mollona e Presutti(2006)

La competizione si struttura nel seguente modo: x1, x2 fornitori dello stesso distrettocompetono tra loro per cercare per cercare di rifornire il maggior numero di produttorirappresentati da y1 +y2. A loro volta i produttori y1, y2 dello stesso distretto competonoper cercare di accaparrarsi il maggior numero di consumatori finali u. Si definisce quindiun sistema omogeneo nella seguente forma:

dx1

dt= g1x1(

(y1 + y2)− c21 − x1

(y1 + y2))

dx2

dt= g2x2(

(y1 + y2)− c12 − x2

(y1 + y2))

dy1

dt= giy1(

u− c21 − y1

u)

dy2

dt= giy2(

u− c12 − y2

u)

49

3.3 Analisi del Fenomeno Empirico: Ricerca dei Pun-

ti Critici

Definito il modello di riferimento cercheremo di studiare qualitativamente la dinamicaevolutiva dello stesso; sia da un punto di vista analitico, cercheremo i punti critici enei loro intorni analizzeremo la stabilita sia da un punto di vista computazionale ossiacercando di capire cosa succede al sistema quando ci si allontana da questi punti.

dx1

dt= g1x1(

(y1 + y2)− c21 − x1

(y1 + y2))

dx2

dt= g2x2(

(y1 + y2)− c12 − x2

(y1 + y2))

dy1

dt= g1y1(

u− c21 − y1

u)

dy2

dt= g2y2(

u− c12 − y2

u)

Dove:

• x1, x2, y1 e y2 sono variabili che rappresentano le popolazioni

• g1, g2, g1 e g2 sono costanti positive che rappresentano i tassi di crescita relativi aogni popolazione

• c21, c12, c21 e c12 sono costanti positive che rappresentano i tassi di competitivita

• u e una costante positiva che rappresenta l’ampiezza del mercato cioe il numerodei potenziali clienti

Si osserva che le ultime due equazioni del sistema formano un sistema indipendentedel tipo “Competizione tra Specie ”; per cui e molto semplice calcolare subito le possibilicoordinate y1, y2 dei punti critici. Sostituendo questi valori nelle prime due equazionericaveremo le restanti coordinate di x1 e x2.

50

Il sistema ammette dodici punti di equilibrio, cioe punti di particolare importanza incui sistema prende comportamenti particolari di stabilita o instabilita.

P1 = (0, 0, 0, u)

P2 = (0, u, 0, u)

P3 = (u, 0, 0, u)

P4 = (u(c21 − 1)

c12c21− 1,u(c12 − 1)

c12c21− 1, 0, u)

P5 = (0, 0, u, 0)

P6 = (0, u, u, 0)

P7 = (u, 0, u, 0)

P8 = (u(c21 − 1)

c12c21 − 1,u(c12 − 1)

c12c21 − 1, u, 0)

P9 = (0, 0,u(1− c21)

1− c12c21,u(1− c12)

1− c12c21)

P10 = (0,u(2− c12 − c21)

1− c12c21

,u(1− c21)

1− c12c21

,u(1− c12)

1− c12c21

)

P11 = (u(2− c12 − c21)

1− c12c21

, 0,u(1− c21)

1− c12c21

,u(1− c12)

1− c12c21

)

P12 = (u(2− c12 − c21)(1− c21)

1− c12c21(1− c12c21),u(2− c12 − c21)(1− c12)

1− c12c21(1− c12c21),u(1− c21)

1− c12c21

,u(1− c12)

1− c12c21

)

51

3.4 Linearizzazione del Fenomeno Empirico

Per caratterizzare i punti di equilibrio e darne una descrizione occorre studiarne la sta-bilita. Il problema in esame e definito da un sistema non lineare autonomo compostoda quattro equazioni differenziali. Per maggiore semplicita chiameremo l’equazione dif-ferenziale i-esima con Fi e la derivata della funzione i-esima nel punto P0 rispetto allavariabile j con F i

j (P0).Per comprendere un sistema cosı complesso applicheremo il metodo indiretto di Lya-punov o linearizzazione, cioe trasleremo nell’origine i nostri punti critici annullando laparte non non lineare.

Calcolando le derivate parziali otteniamo

F 1x1

= g1 −g1c21x2

y1 + y2

− 2g1x1

y1 + y2

F 1x2

= −c12g1x1

y1 + y2

F 1y1

=g1x1(x1 + c21x2)

y1 + y2

F 1y2

=g1x1(x1 + c21x2)

(y1 + y2)2

F 2x1

= −c12g2x2

y1 + y2

F 2x2

= g2 −g2c12x1

y1 + y2

− 2g2x2

y1 + y2

52

F 2y1

=g2x2(x2 + c12x1)

(y1 + y2)2

F 2y2

=g2x2(x2 + c12x1)

(y1 + y2)2

F 3x1

= 0

F 3x2

= 0

F 3y1

=g1(u− 2y1 − c12y2)

u

F 3y2

= − c12g1y1

u

F 4x1

= 0

F 4x2

= 0

F 4y1

= − c12g2y2

u

F 4y2

=g2(u− 2y2 − c12y1)

u

e quindi il sistema lineare approssimazione del sistema non lineare, sara ottenutosostituendo le derivate alla matrice cosı costruita.

F 1x1

(P0) F 1x2

(P0) F 1y1

(P0) F 1y2

(P0)F 2x1

(P0) F 2x2

(P0) F 2y1

(P0) F 2y2

(P0)F 3x1

(P0) F 3x2

(P0) F 3y1

(P0) F 3y2

(P0)F 4x1

(P0) F 4x2

(P0) F 4y1

(P0) F 4y2

(P0)

53

Otteniamo quindi il seguente sistema linearizzato.

d

dt

x1

x2

y1

y2

=

F 1x1

(P0) F 1x2

(P0) F 1y1

(P0) F 1y2

(P0)F 2x1

(P0) F 2x2

(P0) F 2y1

(P0) F 2y2

(P0)0 0 F 3

y1(P0) F 3

y2(P0)

0 0 F 4y1

(P0) F 4y2

(P0)

x1

x2

y1

y2

Si osserva immediatamente, che la matrice dei coefficienti e composta da quattro sot-tomatrici di ordine due, di cui una, quella in basso a sinistra ha tutti elementi nulli.In questo caso particolare e possibile determinare gli autovalori o strade del polinomiocaratteristico det(A − Ir) della matrice principale semplicemente risolvendo le due sot-tomatrici di ordine secondo sulla diagonale che chiameremo D1 e D2 in un generico puntoP0.

D1 =

(F 1x1

(P0) F 1x2

(P0)F 2x1

(P0) F 2x2

(P0)

)=

(g1 − g1c21x0

2

y01+y02− 2g1x0

1

y01+y02− c21g1x0

1

y01+y02

− c12g2x02

y01+y02g2 − g2c12x0

1

y01+y02− 2g2x0

2

y01+y02

)

D2 =

(F 3y1

(P0) F 3y2

(P0)F 4y1

(P0) F 4y2

(P0)

)=

(g1(u−2y01−c12y02)

u− c21g1y01

u

− c12g2y02u

g2(u−2y02−c12y01)

u

)

Si osserva che ciascun blocco si comporta come un modello “Competizione tra Specie”, ed e utile notare che la stabilita del blocco D1 dipende direttamente dal livello inferioreossia dal blocco D2 che invece agisce in totale indipendenza. Grazie a queste deduzionivedremo che per calcolare la stabilita del sistema non sara piu necessario calcolare gliautovalori ma bastera applicare le formule del modello “Competizione tra Specie ”giaampiamente discusse.

ε1 = g1 , ε2 = g2 , σ1 = g1/(y01 + y0

2) , σ2 = g2/(y01 + y0

2)α1 = c21g1/(y

01 + y0

2) , α2 = c12g2/(y01 + y0

2) ,ε1 = g1 , ε2 = g2 ,

σ1 = g1/u , σ2 = g2/u ,α1 = c21g1/u , α2 = c12g2/u ,

54

Riprendendo la teoria delle “Competizione tra Specie ”, vediamo che un punto criticoe stabile se σ1σ2 − α1α2 > 0, Quindi la coesistenza sara possibile solo in condizioni dilimitata interazione. Sostituendo alla condizione il livello dei fornitori otteniamo

gx1

(y01 + y0

2)

gx2

(y01 + y0

2)− cx2x1gx1

(y01 + y0

2)

cx1x2gx2

(y01 + y0

2)> 0

gx1gx2

(y01 + y0

2)2(1− cx2x1cx1x2) > 0

da cui P0 e stabile se c12c21 < 1; allo stesso modo per il livello dei produttori otteniamoc21c12 < 1. La stabilita del sistema e quindi definita dalla stabilita delle due sotto matriciD1, D2 e dalle loro relative condizioni. Se solo una o entrambe non sono verificate allorail punto di equilibrio e instabile.

3.5 Analisi Qualitativa: Fenomeno Empirico

Dopo aver osservato le condizioni applicate al nostro sistema analizziamo la qualita diogni punto critico. Sapendo che (0,0) e sempre instabile nel modello base, i punti criticiP1, P5, P9 saranno a loro volta sempre instabili. Osserviamo che nella sotto matrice D1

det(D1 − Ir) = det(

(F 1x1

(P0) F 1x2

(P0)F 2x1

(P0) F 2x2

(P0)

)− Ir) = det(

(g1 − r 0

0 g2 − r

))

gli autovalori saranno sempre positivi portando instabilita a tutto il sistema nel punto.

55

Proseguiamo nell’analisi per i punti critici P2, P3, P4, P7, P8, P10, P11, P12, partendo dalcentro del sistema, cioe dal punto P12 e successivamente tipicizzando tutti gli altri.

Riprendendo i concetti del modello “Competizione tra Specie ”diciamo che:

• Se il punto con coordinate positive e stabile, i punti che si trovano sugli assi, cioequei punti con una coordinata nulla, sono instabili.

• Se il punto con coordinate positive e instabile, similarmente al caso precedente ipunti sugli assi sono stabili.

Supponendo P12 sia stabile e quindi c12c21 < 1 e c12c21 < 1 verificate, identificheremola natura degli altri punti dalle condizioni ε1σ2 − α1ε2 > 0 e σ1ε2 − α2ε1 > 0.

Sostituendo otterremo le seguenti condizioni

g1g2

(y01 + y0

2)− g1g2c21

(y01 + y0

2)> 0 → g1g2

(y01 + y0

2)(1− c21) > 0 → c21 < 1

g1g2

(y01 + y0

2)− g1g2c12

(y01 + y0

2)> 0 → g1g2

(y01 + y0

2)(1− c12) > 0 → c12 < 1

g1g2

u− g1g2c21

u> 0 → g1g2

u(1− c21) > 0 → c21 < 1

g1g2

u− g1g2c12

u> 0 → g1g2

u(1− c12) > 0 → c12 < 1

dalle quali osserviamo che tutti i punti con una coordinata nulla per le osservazioniprecedenti sono instabili.

56

Analogamente per un punto P12 instabile, riapplicheremo lo stesso procedimentoidentificando tre casistiche per gli altri punti.

• c12c21 > 1 e c21c21 < 1: l’instabilita del punto P12 e data dalle ultime due equazionicioe dalle ultime due coordinate del punto, mentre le prime due sono stabili. Perla positivita delle coordinate del punto P12, vale per gli altri punti che:

c12 > 1 e c21 > 1, c21 < 1 e c21 < 1, cioe i punti P4, P8 sono stabili per effettodella coordinata nulla, mentre i restanti P2, P3, P6, P7, P10, P11 no.

• c12c21 > 1 e c21c21 > 1: l’instabilita del punto P12 e data da tutte le equazioni cioeda entrambe le coppie di coordinate. Per la positivita delle coordinate del puntoP12, vale per gli altri punti che:

c12 > 1 e c21 > 1, c21 > 1 e c21 > 1, cioe risultano stabili i pun-ti P2, P3, P6, P7, P10, P11 per effetto della coordinata nulla in almeno una coppia,mentre P4, P8 risultano instabili.

• c21c21 > 1 e c12c21 < 1: l’instabilita del punto P12 e determinata dalle prime dueequazioni e conseguentemente dalle prime due coordinate. Per la positivita dellecoordinate del punto P12, vale per gli altri punti che:

c21 > 1 e c21 > 1, c12 > 1 e c21 < 1, cioe risultano stabili i punti P10, P11 pereffetto della coordinata nulla, mentre P2, P3, P4, P6, P7, P8 risultano instabili.

57

3.6 Significato Economico dei Punti Critici

Dato un punto di equilibrio generico P0 = (x1, x2, y1, y2), cercheremo ora di identifi-care tutti quegli aspetti economici che l’analisi matematica non ha messo in evidenza.Riprendendo i concetti precedentemente visti, si ricorda che ogni coordinata di un puntogenerico corrisponde e rappresenta la numerosita di una delle popolazioni di fornitori inesame nel nostro distretto. Come gia detto le popolazioni xi competono tra loro cosıcome le popolazioni yi.

Inizieremo il nostro studio osservando i punti definiti dall’analisi matematica comesempre instabili; cioe P1, P5, P6, P9. Tutti questi punti sono accomunati da delle coor-dinate xi nulle, responsabili dell’instabilita del sistema nel punto. Da una prospettivaeconomica esse comportano la mancanza dei fornitori. Fortunatamente possiamo pensarea una fornitura di prodotti esterna al sistema che altrimenti sarebbe privo di significa-to. Per tutti gli altri punti, non sara possibile in generale darne una tipicizzazione. Lospazio di appartenenza del nostro sistema non e rappresentabile e pertanto di difficilecomprensione per quanto concerne le traiettorie. Osserviamo ad ogni modo che scostan-doci dall’origine, anche di poco, il sistema iniziera una evoluzione, regolata e definitadai tassi di competitivita. Ogni punto avra una natura stabile o instabile al verificarsi dicerte condizioni di competitivita fra le popolazioni. Grazie all’analisi svolta, si sono iden-tificati alcuni gruppi di appartenenze per i punti sotto certe condizioni, cioe ogni puntonella stessa condizione economica rispondera allo stesso modo. Considerando i puntiP2, P3, P7, caratterizzati dalla presenza di una sola popolazione di fornitori e da una soladi produttori, i suddetti punti possono essere contemporaneamente stabili e quindi strut-turare diverse situazioni economiche tali per cui non verranno accettati nuovi entranti.D’altra parte, in un contesto diametralmente opposto, essi potranno essere instabili equindi incentivare l’ingresso di nuovi competitori che sposteranno il sistema verso nuovipunti di equilibrio. I punti P4, P8, sono caratterizzati dalla presenza di una sola popo-lazione di produttori, ne consegue che in mancanza dell’altra, essa dominera il mercatodei potenziali clienti. E’ importante aggiungere che in un punto di equilibrio con unacoordinata nulla, potenzialmente a seconda della competitivita, potranno emergere o nonemergere nuovi competitori stranieri. Invece nei punti P10, P11, avremo una situazionecontraria rispetto al caso precedente; cioe la presenza di una sola popolazione di forni-tori che potra rifornire entrambi i produttori. Infine, il punto P12, e caratterizzato dallapresenza positiva di tutte le popolazioni, che a prescindere dal numero, non potrannoalterarne le caratteristiche di stabilita. La spiegazione del fenomeno e dovuta all’unicofattore rilevante per il sistema: la competitivita delle diverse popolazioni.

58

3.7 Funzione di Lyapunov e Stabilita del Fenomeno

Empirico

Terminata l’analisi qualitativa nelle due sezioni precedenti, osserviamo la totale mancan-za di informazioni riguardanti i bacini di attrazioni che invece rivestono un importanzafondamentale in quanto ci dicono di quanto possiamo allontanarci da un punto stabilesenza che il sistema non vi ritorni. Ricordiamo che un sistema e stabile se la sua dinam-ica, partendo da certe condizioni iniziali in un certo intorno del punto, raggiunge talepunto. Puo essere interessante comprendere gli spostamenti massimali possibili, cioe diquanto possono variare le costanti senza allontanarsi dal punto. Analizzando il problemada un differente punto di vista, si cerchera di identificare prima una funzione di Lyapunovper il centro del sistema cioe il P12 e successivamente un dominio adeguato.

L’applicazione del metodo diretto di Lyapunov comporta svariate difficolta, tra questala maggiore e sicuramente l’assenza di un metodo generico per creare la funzione scalareenergia necessaria all’analisi della stabilita del sistema.Ricordando le nozioni introdotte nel secondo capitolo sulla costruzione di una funzioneenergia, dobbiamo trovare una funzione V , definita positiva in un qualche dominio Dtale che la sua derivata V sia definita negativa nello stesso dominio. Richiamando ilmetodo di costruzione utile per funzioni definite positive.

V (x, y) = ax2 + bxy + cy2

• e positiva definita se e solo se a > 0 e 4ac− b2 > 0

• e negativa definita se e solo se a < 0 e 4ac− b2 < 0

Osserviamo che ponendo b = 0, a > 0 e c > 0 la condizione di positivita definita esempre rispettata e la funzione e nella forma.

V (x, y) = ax2 + cy2

59

Utilizzando come base l’esempio del capitolo due e il metodo appena dimostrato,ottenniamo una funzione di Lyapunov candidata in quattro variabili del tipo:

V (α, β, γ, δ) = 1/2α2 + 1/2β2 + 1/2γ2 + 1/2δ2

ove a = 1/2, b = 1/2, c = 1/2 e d = 1/2 tutte positive, scelte per semplificare lederivate della funzione. Osserviamo che chiaramente V e definita positiva in D, cioeV (x1, x2, y1, y2) > 0 per tutti gli altri punti in D e V (0, 0, 0, 0) = 0. Supponendo unpunto critico generico xeq = (α, β, γ, δ), isolato nell’origine e definito sostituendo alleequazioni differenziali. Da cui:

• x1 = X1 + α

• x2 = X2 + β

• y1 = Y1 + γ

• y2 = Y2 + δ

60

e quindi utilizzando il programma Mathematica per svolgere i calcoli simbolicamenteotteniamo:

x1=(u(2−cc12−cc21)(1−c21))/((1−cc12∗cc21)(1−c12∗c21))+α

(1−c21)(2−cc12−cc21)u(1−c12c21)(1−cc12cc21)

+α

x2=(u(2−cc12−cc21)(1−c12))/((1−cc12∗cc21)(1−c12∗c21))+β

(1−c12)(2−cc12−cc21)u(1−c12c21)(1−cc12cc21)

+β

y1=(u(1−cc21))/(1−cc12∗cc21)+γ

(1−cc21)u1−cc12cc21+γ

y2=(u(1−cc12)/(1−cc12∗cc21))+δ

(1−cc12)u1−cc12cc21+δ

ove (X1, X2, Y1, Y2), sono le coordinate generiche del punto critico. La sua derivataparziale nel punto P12 centro del sistema, sara definita nel seguente modo:

V (α, β, γ, δ) =

Vα(α, β, γ, δ)dx1

dt+ Vβ(α, β, γ, δ)

dx2

dt+ Vγ(α, β, γ, δ)

dy1

dt+ Vδ(α, β, γ, δ)

dy2

dt

da cui otteniamo

61

VV=α∗dX1+β∗dX2+γ∗dY1+δ∗dY2

gg2(u− (1−cc12)u1−cc12cc21−cc12( (1−cc21)u

1−cc12cc21+γ)−δ)δ( (1−cc12)u1−cc12cc21+δ)

u +

g2β( (1−c12)(2−cc12−cc21)u(1−c12c21)(1−cc12cc21)

+β)( (1−cc12)u1−cc12cc21+

(1−cc21)u1−cc12cc21−

(1−c12)(2−cc12−cc21)u(1−c12c21)(1−cc12cc21)

−c12( (1−c21)(2−cc12−cc21)u(1−c12c21)(1−cc12cc21)

+α)−β+γ+δ)(1−cc12)u1−cc12cc21+

(1−cc21)u1−cc12cc21+γ+δ

+

g2β( (1−c12)(2−cc12−cc21)u(1−c12c21)(1−cc12cc21)

+β)( (1−cc12)u1−cc12cc21+

(1−cc21)u1−cc12cc21−

(1−c12)(2−cc12−cc21)u(1−c12c21)(1−cc12cc21)

−c12( (1−c21)(2−cc12−cc21)u(1−c12c21)(1−cc12cc21)

+α)−β+γ+δ)(1−cc12)u1−cc12cc21+

(1−cc21)u1−cc12cc21+γ+δ

+

g1α( (1−c21)(2−cc12−cc21)u(1−c12c21)(1−cc12cc21)

+α)( (1−cc12)u1−cc12cc21+

(1−cc21)u1−cc12cc21−

(1−c21)(2−cc12−cc21)u(1−c12c21)(1−cc12cc21)

−α−c21( (1−c12)(2−cc12−cc21)u(1−c12c21)(1−cc12cc21)

+β)+γ+δ)(1−cc12)u1−cc12cc21+

(1−cc21)u1−cc12cc21+γ+δ

+

gg1γ( (1−cc21)u1−cc12cc21+γ)(u− (1−cc21)u

1−cc12cc21−γ−cc21( (1−cc12)u1−cc12cc21+δ))

u

Per maggiore chiarezza, ricordiamo il primo teorema di Lyapunov:

“Si suppone che il sistema autonomo abbia un punto critico isolato all’origine. Seesiste una funzione V continua detta funzione energia con derivate parziali del primoordine continue positivamente definite e per le quali V (x, y) definita comeVx(x, y)dx

dt+ Vy(x, y)dy

dte negativamente definita in un qualche dominio D nel piano-xy

contenente (0,0), allora l’origine e asintoticamente stabile.Se V e negativa ma solamente semi-definita allora l’origine e un punto critico stabile.”

Il passo successivo e la dimostrazione che la nostra funzione candidata V , e unafunzione di Lyapunov, cioe soddisfa le ipotesi di stabilita introdotte dal primo teore-ma. Per rispondere a questo quesito, determineremo per prima cosa il gradiente di V esuccessivamente l’hessiano, entrambi nell’origine del nostro piano αβγδ.

62

Calcoliamo il gradiente nell’origine per il punto critico P12. Ricordiamo che se f :X → R, con X aperto di Rn generica in n variabili nel punto (x1, x2, ..., xn) e un vettoredi Rn avente per componenti le derivate parziali prime calcolate nel punto,

~∇ f(α0, β0, γ0, δ0) =

fα(α0, β0, γ0, δ0)fβ(α0, β0, γ0, δ0)fγ(α0, β0, γ0, δ0)fδ(α0, β0, γ0, δ0)

Nel nostro caso possiamo calcolarle usando il programma Mathematica:

D[VV,α]/.{α→0,β→0,γ→0,δ→0}

(1−c21)(2−cc12−cc21)g1u( (1−cc12)u1−cc12cc21+

(1−cc21)u1−cc12cc21−

(1−c21)(2−cc12−cc21)u(1−c12c21)(1−cc12cc21)

− (1−c12)c21(2−cc12−cc21)u(1−c12c21)(1−cc12cc21) )

(1−c12c21)(1−cc12cc21)( (1−cc12)u1−cc12cc21+

(1−cc21)u1−cc12cc21)

Simplify[%]

0

D[VV,β]/.{α→0,β→0,γ→0,δ→0}

(1−c12)(2−cc12−cc21)g2u( (1−cc12)u1−cc12cc21+

(1−cc21)u1−cc12cc21−

(1−c12)(2−cc12−cc21)u(1−c12c21)(1−cc12cc21)

− c12(1−c21)(2−cc12−cc21)u(1−c12c21)(1−cc12cc21) )

(1−c12c21)(1−cc12cc21)( (1−cc12)u1−cc12cc21+

(1−cc21)u1−cc12cc21)

Simplify[%]

0

D[VV,γ]/.{α→0,β→0,γ→0,δ→0}

(1−cc21)gg1(u− (1−cc21)u1−cc12cc21−

(1−cc12)cc21u1−cc12cc21 )

1−cc12cc21

Simplify[%]

0

D[VV,δ]/.{α→0,β→0,γ→0,δ→0}

(1−cc12)gg2(u− (1−cc12)u1−cc12cc21−

cc12(1−cc21)u1−cc12cc21 )

1−cc12cc21

Simplify[%]

0

ove osserviamo che ~∇V = 0 e quindi in quel punto la funzione e stazionaria. Verificataquesta prima condizione necessaria, possiamo procedere alla dimostrazione dell’ipotesidi stabilita.

63

Calcoliamo la matrice hessiana della funzione V nel punto critico, definita per una

funzione generica di n variabili f(x) = f(x1, x2, ..., xn) come Hf (x)ij = ∂2f(x)∂xi∂xj

, da cui:

HV (α0, β0, γ0, δ0) =

∂ f2

∂α2∂ f2

∂αβ∂ f2

∂αγ∂ f2

∂αδ∂ f2

∂βα∂ f2

∂β2∂ f2

∂βγ∂ f2

∂βδ∂ f2

∂γα∂ f2

∂γβ∂ f2

∂γ2∂ f2

∂γδ∂ f2

∂δα∂ f2

∂δβ∂ f2

∂δγ∂ f2

∂δ2

e quindi nel nostro caso:

MatrixForm[HH]

− 2(−1+c21)g1−1+c12c21 −r

c21g1−c212g1−(−1+c12)c12g2−1+c12c21

(−1+c21)g1−1+c12c21

(−1+c21)g1−1+c12c21

c21g1−c212g1−(−1+c12)c12g2−1+c12c21 − 2(−1+c12)g2

−1+c12c21 −r(−1+c12)g2−1+c12c21

(−1+c12)g2−1+c12c21

(−1+c21)g1−1+c12c21

(−1+c12)g2−1+c12c21 − 2(−1+cc21)gg1

−1+cc12cc21 −rcc21gg1−cc212gg1−(−1+cc12)cc12gg2

−1+cc12cc21

(−1+c21)g1−1+c12c21

(−1+c12)g2−1+c12c21

cc21gg1−cc212gg1−(−1+cc12)cc12gg2−1+cc12cc21 − 2(−1+cc12)gg2

−1+cc12cc21 −r

64

Calcoliamo il discriminante con Mathematica:

detHH=Simplify[Det[HH]]

detHH=Simplify[Det[HH]]

((2−3cc21+cc212)gg1+cc122gg2+r−cc12(gg2+cc21r))+

1−1+cc12cc21 (−1+c12)g2((2−3cc21+cc212)gg1+cc122gg2+r−cc12(gg2+cc21r))

(2c212g12−c213g12+(−1+c12)g2((−2+c12)g1−r)−c21(g12+(2−3c12+c122)g1g2−(−1+c12)c12g2r))−

1−1+cc12cc21(cc21gg1−cc212gg1−(−1+cc12)cc12gg2)

(−(−1+c21)g1((−1+c12)2c12g22+(−1+c21)g1((−1+c12)(−2+c21)g2+r−c12c21r))+

1−1+cc12cc21 (−1+c12c21)(cc21gg1−cc212gg1−(−1+cc12)cc12gg2)

(−(−c21g1+c212g1+(−1+c12)c12g2)

2+

(2(−1+c21)g1+(−1+c12c21)r)(2(−1+c12)g2+(−1+c12c21)r))−(−1+c12)g2

(−2c212g12+c213g12−(−1+c12)g2((−2+c12)g1−r)+c21(g12+(2−3c12+c122)g1g2−(−1+c12)c12g2r)))+

(− 2(−1+cc12)gg2−1+cc12cc21 −r)(−(−1+c21)g1((−1+c12)2c12g22+(−1+c21)g1((−1+c12)(−2+c21)g2+r−c12c21r))+

(−1+c12c21)(− 2(−1+cc21)gg1−1+cc12cc21 −r)

(−(−c21g1+c212g1+(−1+c12)c12g2)

2+

(2(−1+c21)g1+(−1+c12c21)r)(2(−1+c12)g2+(−1+c12c21)r))−(−1+c12)g2

(−2c212g12+c213g12−(−1+c12)g2((−2+c12)g1−r)+c21(g12+(2−3c12+c122)g1g2−(−1+c12)c12g2r))))

Vediamo che dato il numero di variabili, non e possibile definire una legge che regolila negativita degli autovalori. Di conseguenza per valutare se la funzione V e definitanegativa, sostituiamo valori numerici alle costanti.

• u=1000

• gx1=gx2=gy1=gy2=0,05

• cx2x1=0.6, cx1x2=0.2, cy2y1=0.8, cy1y2=0.4

Osservando immediatamente che la condizione di stabilita per il punto P12 σ1σ2 −α1α2 > 0 e rispettata, conseguentemente la funzione V deve avere un minimo localenell’origine.

Ricalcolando il discriminante dopo la sostituzione otteniamo:

,

65

Figura 3.2: Risultati Matrice Hessiana in P12

quattro autovalori negativi, cioe la funzione nell’origine ha un massimo. Questoimportante risultato dimostra che effettivamente V esiste ed e una funzione di Lyapunove che sicuramente in un certo intorno D dell’origine, essa e negativa, quindi il puntocritico P12 e stabile. Sfortunatamente con il metodo analitico e data la complessitadel fenomeno in esame non siamo riusciti a determinare un dominio per il bacino diattrazione, che ad ogni modo puo essere determinato attraverso il calcolo numerico.

66

Bibliografia

[1] William E., Boyce Richard C. DiPrima: Elementary Differential Equations andBoundary Value Problems, (2001). Ed. Wiley.

[2] Mollona E., M. Presutti: Proceedings of the Twenty-fourth International Sys-tem Dynamics Conference, (2006). ALBANY: System Dynamics Society (UNITEDSTATES).

[3] Tesi Paparoni

67

Ringraziamenti