#!/usr/bin/python
from pylab import *
#r = 3.83
#x0 = 0.1
#N = 100
def orbita(r,x0,N):
"""
Retorna un vector de dimensio N amb els N primers iterats per aplicacio amb seed x0 i parametre r
x=orbita(3.8,0.1,100)
"""
x = [x0]
for n in xrange(0,N):
x.append( aplicacio(r,x[n]) )
return x
def dibuixaorbita(r,x0,N):
"""
Retorna el dibuix amb els N primers iterats per aplicacio amb seed x0 i parametre r.
Exemple: dibuixaorbita(3.8,0.1,100); show()
"""
xlabel('n') # set x-axis label
ylabel('x(n)') # set y-axis label
title('r= ' + str(r)) # set plot title
x=orbita(r,x0,N)
p=plot(x, 'ro', x , 'b')
return p
def duesorbites(r,x0,N):
"""
Retorna dos vectors de dimensio N amb els N primers iterats per aplicacio amb seed x0 (vector de 2 comp) i parametre r
Exemple:
x1,x2=duesorbites(3.8,[0.1,0.1+1e-10],1000)
"""
x1 = [x0[0]]
x2 = [x0[1]]
for n in xrange(0,N):
x1.append( aplicacio(r,x1[n]) )
x2.append( aplicacio(r,x2[n]) )
return x1,x2
def dibuixaduesorbites(r,x0,N):
"""
Retorna el dibuix de dues orbites amb els N primers iterats per aplicacio amb seed x0 (vector de 2 comp) i parametre r. Es util per veure la sensitivitat de ci.
dibuixaduesorbites(3.8,[0.1,0.1+1e-10],100); show()
"""
x1,x2=duesorbites(r,x0,N)
xlabel('n')
ylabel('x(n)')
title(' r= ' + str(r))
p=plot(x1, 'ro', x1 , 'b')
p+=plot(x2, 'yD', x2 , 'g')
return p
def aplicacio(r,x):
"""
Retorna laplicacio logistica:
"""
return r * x * (1.0 - x)
def dibuixaplicacio(r,iteracions):
"""
Retorna el dibuix dels iterats de laplicacio amb parametre r
exemple: dibuixaplicacio(3.8,4); show()
"""
x = linspace(0.0,1.0,1000)
TitleString = 'r=%g' % r
title(TitleString)
xlabel('x(n)')
ylabel('x(n+1)')
p=plot(x, x, 'b--')
#logistica
y=x
for n in arange(0,iteracions):
y=aplicacio(r,y)
p+=plot(x, y)
return p
def cobweb(r,inicial,iteracions):
"""
Retorna dos vectors x0 i x1 de dimensio iteracions que contenen les x i les y del diagrama de teranyina de laplicacio.
exemple: x0,x1=cobweb(3.8, 0.1, 10)
"""
state = inicial
x0 = [ ] # Valor de x_n
x1 = [ ] # Valor de x_{n+1}
# Anem actualitzant els vectors
for n in xrange(iteracions):
x0.append( state )
x1.append( state )
x0.append( state )
state = aplicacio(r,state)
x1.append( state )
return x0,x1
def dibuixacobweb(r,inicial,iteracions):
"""
Retorna el dibuix de la teranyina amb parametre r, condicions incials inicial, i iteracions.
exemple: dibuixacobweb(3.8, 0.1, 10); show()
o be: dibuixacobweb(3.8, 0.1, 10)+dibuixaplicacio(3.8,1); show()
"""
x0,x1= cobweb(r,inicial,iteracions)
p=plot(x0, x1, '+-')
return p
def dibuixabifurcacio(rlow,rhigh,iteratsx = 1100,transientx = 1000,samplesr = 1000):
"""
Dibuixa el diagrama de bifurcacio de l'aplicacio entre rlow i rhigh. Els arguments opcionals son:
- iteratsx = 1100, nombre d'iterats de cada orbita
- transientx = 1000, transient que deixem passar abans de dibuixar
- samplesr = 1000, mostres en el parametre r.
exemple: dibuixabifurcacio(2.56,4)
o be: dibuixabifurcacio(2.56,4,110,100,10000)
"""
figure(1,(8,6))
TitleString = 'r a [%g,%g]' % (rlow,rhigh)
title(TitleString)
xlabel('r')
ylabel('x(n)')
# Diem al plot() que autoscali en el rang que li marquem per rhigh i rlow
plot([rhigh], [1.0], 'k,')
plot([rhigh], [0.0], 'k,')
plot([rlow], [0.0], 'k,')
plot([rlow], [1.0], 'k,')
#posem una condicio inicial
ic = 0.2
#Definim l'increment per a la x
rInc = (rhigh-rlow)/float(samplesr)
for r in arange(rlow,rhigh,rInc):
# Definim la ci
state = ic
# En el transientx no registrem res
for i in xrange(transientx):
state = aplicacio(r,state)
# A partir d'ara ja gravem els valors
rsweep = [ ] # valor del parametre
x = [ ] # iterats
for i in xrange(iteratsx):
state = aplicacio(r,state)
rsweep.append(r)
x.append( state )
plot(rsweep, x, 'b,') # Dibuixem com a pixels vermells
show()4 L’aplicació logística
4.1 Estudi de les bifurcacions en el model logístic
El fitxer logistica.py conté les definicions per estudiar les bifurcacions de l’aplicació logística. Per obtenir la informació de les funcions podem fer, des de la consola d’una terminal
pydoc -w ./logistica.py
ens crearà un fitxer anomenat logistica.html amb tota la documentació, no només de les funcions que hem creat sinó també de les altres. El contingut del fitxer és el següent:
Reproduiu l’anàlisi amb l’aplicació de Ricker: \[x_{n+1}= a exp(-x_n), \quad n \ge 0\] que no té el problema que doni valors negatius de la població. Veieu que, al començament, té únic punt fix estable que bifurca a \(a=8\) a una òrbita periòdica de periode 2. Al seu torn aquest bifurca a \(a=12\). Construïu la cascada de bifurcacions, tal com hem fet amb la logística.
4.2 Càlcul dels exponents de Lyapunov
El fitxer lyapunov.py conté el càlcul i dibuix dels exponents de Lyapunov:
from scipy import arange,log
from pylab import *
from logistica import *
# Aqui importem totes les funcions
# definim l'aplicacio derivada
def aplicacioderivada(r,x):
return r*(1.0-2.0*x)
#valors del rastrejat
plow = 0.1
phigh = 4.0
# Aquest sera el numero de passos de la a
nSteps = 1000.0
ic=0.2
# Transient dels iterats (abans de calcular)
nTransients = 100
# Nombre total t'iteracions
nIterates = 2500
# Condicio inicial de l'orbita
ic = 0.2
# Comencem a dibuixar
TitleString = 'Exponent de Lyapunov per a [%g,%g]' % (plow,phigh)
title(TitleString)
# Eixos
xlabel('r')
ylabel('$\lambda$')
# \lambda=0 representa la transicio
axhline(0.0,0.0,1.0,color = 'k')
# A psweep guardarem els valors de r i a lce el valor de l'exponent de Lyapunov
psweep = arange(plow,phigh,(phigh-plow)/nSteps)# Parametres
lce = [ ] # Exponents de Lyapunov
# Anem escanejant el
for p in psweep:
# definim la condicio inicial
x = ic
# Ens oblidem de les primeres iteracions
for i in xrange(nTransients):
x = aplicacio(p,x)
# Anem a estimar l'exponent de Lyapunov amb els restants
# incialitzem a zero l'exponent de Lyapunov
l = 0.0
for i in xrange(nIterates):
x = aplicacio(p,x)
l = l + log(abs(aplicacioderivada(p,x)))
# dividim pel nombre d'iterats
l = l / float(nIterates)
lce.append(l)
# Ho dibuixem tot
plot(psweep, lce, 'r-')
# Ho mostrem
show()Observem com crida al fitxer logistica.py.
Calculeu els exponents de Lyapunov per a l’aplicació de Ricker i compareu-los amb el diagrama de bifurcació que hem obtingut a l’altre apartat. En quines àrees sembla haver-hi caos?
4.3 Un model per a la vespa
Barlow i altres (bio:wasp-ricker?), l’any 1992 van estudiar la vespa comuna (vespula vulgaris) que s’havia introduït a Nova Zelanda a partir dels anys 1980. Aquesta espècia és força agressiva amb els altres insectes autòctons i amb els animals. Els autors van registrar les poblacions anuals a la primavera (en termes de nius per hectàrea) en 7 llocs diferents del 1988 al 1992, tal i com s’expressa a la taula següent.
| Lloc | 1988 | 1989 | 1990 | 1991 | 1992 |
| 1 | 8.6 | 31.1 | 7.0 | 11.7 | 10.2 |
| 2 | 2.7 | 6.9 | 3.3 | 4.4 | 3.1 |
| 3 | 10.5 | 15.8 | 8.2 | 11.6 | 12.1 |
| 4 | 1.0 | 1.9 | 6.0 | 1.0 | 1.0 |
| 5 | 16.0 | 11.8 | 10 | 15.7 | 19.9 |
| 6 | 18.5 | 32.9 | 17.1 | 13.6 | 13.0 |
| 7 | 20.5 | 19.8 | 12.9 | 15.7 | 10.6 |
Preneu les dades del lloc 1 i intenteu trobar els paràmetres \(a\) i \(b\) que millor s’ajusten a suposar que \[x_{n+1}= a x_n \exp(-bx_n)\] és a dir, que la població ve donada per un model de Ricker. Quin és el comportament qualitatiu?
Com us ho faríeu per estudiar si el comportament (paràmetres \(a\) i \(b\)) són els mateixos per a tots 7 llocs? Podeu fer algun tipus d’anàlisi estadística?