#bloc1.py
from numpy import *
import pylab as p
# Definim els parametres
a = 1.
b = 0.5
c = 1.5
def dX_dt(X, t=0):
""" Retorna les equacions. El t=0 es perque el sistema es autonom """
return array([ X[0]*(1-X[0]-a*X[1]),c*X[1]*(1-b*X[0]-X[1])])5 Models Predador Presa
5.1 Els models de competició
Considerem el model adimensionalitzat de competició: \[\begin{array}{rcl} x' & = & x \left( 1 -x -a y \right) \\ y' & = & c y \left( 1 -b x -y \right) \end{array}\] En primer lloc definim les equacions, tal i com les mostrem en el següent bloc
Ara, des de la línia de comanda, ja podem integrar l’equació:
# bloc2.py
from scipy import integrate
t = linspace(0, 15, 1000) # temps integracio
X0 = array([0.1, 0.8]) # condicions inicials: x=0.1, y=0.8
X, infodict = integrate.odeint(dX_dt, X0, t, full_output=True)
print infodict['message'] # >>> Successful ens dira que ha estat existosPer dibuixar les solucions podem usar el matplotlib:
#bloc3.py
especie1, especie2 = X.T
#Aixo seria el mateix que si fessim:
#especie1=X[:,0]
#especie2=X[:,1]
f1 = p.figure() #iniciem una figura
p.plot(t, especie1, 'r-', label='x (x(0)=%g)' %(especie1[0]))
p.plot(t, especie2, 'b-', label='y (y(0)=%g)' %(especie2[0]))
p.grid() # Ens mostra una graella
p.legend(loc='best') # Ens mostra una llegenda al millor lloc
p.xlabel('$t$')
p.ylabel('poblacio')
p.title("a=%g, b=%g, c= %g" % (a,b,c))
p.show()Observem com podem fer que un escenari i l’altre passin en funció de com siguin \(a\), \(b\), \(c\) i les condicions inicials
5.2 Dibuix del retrat de fase
Com que es tracta d’un sistema de dues dimensions podem intentar comprendre’l millor mitjançant el camp de direccions.
#bloc4.py
#-------------------------------------------------------
# Definim una xarxa i calculem el camp a la xarxa
ymax = 1 # limits dels eixos
xmax = 1
num_punts = 30
x = linspace(0, xmax, num_punts)
y = linspace(0, ymax, num_punts)
X1 , Y1 = meshgrid(x, y) # Creem la xarxa
DX1, DY1 = dX_dt([X1, Y1]) # Taxa de creixement a la xarxa
M = (hypot(DX1, DY1)) # Norma del camp
M[ M == 0] = 1. # Ens evitem que hi hagi divisions per zero
DX1 /= M # Normalitzem
DY1 /= M
#-------------------------------------------------------
# Dibuixem el camp de direccions
# Dibuixem la direccio normalitzada i amb el color de les fletxetes
# la magnitud del camp
p.title("a=%g, b=%g, c= %g" % (a,b,c))
Q = p.quiver(X1, Y1, DX1, DY1, M, pivot='mid', cmap=p.cm.jet)
p.xlabel('x')
p.ylabel('y')
p.legend()
p.grid()
p.xlim(0, xmax)
p.ylim(0, ymax)Si ara volem posar-hi a més unes quantes trajectòries podem fer:
#bloc5.py
valors = linspace(0., 1.0, 20) # variarem les ci entre 0.3 i 0.9
vcolors = p.cm.autumn_r(linspace(0.3, 1., len(valors))) # Per a cada ci triem un color diferent
t = linspace(0, 1000, 10000)
f2 = p.figure()
a=0.25
b=0.5
c=1.5
Xinicial=array([0.0,0.0])
#-------------------------------------------------------
# dibuixem les trajectories que surten dels eixos x
for v, col in zip(valors, vcolors): # Aixo ens fa una doble iteracio
X0 = array([0.01,v]) # condicio inicial
X = integrate.odeint( dX_dt, X0, t) # no posem l'infodict aqui
p.plot( X[:,0], X[:,1], color=col)
# p.plot( X[:,0], X[:,1], lw=3.5*v, color=col, label='X0=(%.f, %.f)' % ( X0[0], X0[1]) ) Per escriure les ci
# dibuixem les trajectories que surten dels eixos y
for v, col in zip(valors, vcolors): # Aixo ens fa una doble iteracio
X0 = array([v,0.01]) # condicio inicial
X = integrate.odeint( dX_dt, X0, t) # no posem l'infodict aqui
p.plot( X[:,0], X[:,1], color=col)
# dibuixem les trajectories que surten dels eixos x
for v, col in zip(valors, vcolors): # Aixo ens fa una doble iteracio
X0 = array([v,0.99]) # condicio inicial
X = integrate.odeint( dX_dt, X0, t) # no posem l'infodict aqui
p.plot( X[:,0], X[:,1], color=col)
# dibuixem les trajectories que surten dels eixos x
for v, col in zip(valors, vcolors): # Aixo ens fa una doble iteracio
X0 = array([0.99,v]) # condicio inicial
X = integrate.odeint( dX_dt, X0, t) # no posem l'infodict aqui
p.plot( X[:,0], X[:,1], color=col)
# dibuixem les trajectories que surten dels eixos x
for v, col in zip(valors, vcolors): # Aixo ens fa una doble iteracio
X0 = array([v,v]) # condicio inicial
X = integrate.odeint( dX_dt, X0, t) # no posem l'infodict aqui
p.plot( X[:,0], X[:,1], color=col)
#-------------------------------------------------------
# Definim una xarxa i calculem el camp a la xarxa
ymax = 1 # limits dels eixos
xmax = 1
num_punts = 30
x = linspace(0, xmax, num_punts)
y = linspace(0, ymax, num_punts)
X1 , Y1 = meshgrid(x, y) # Creem la xarxa
DX1, DY1 = dX_dt([X1, Y1]) # Taxa de creixement a la xarxa
M = (hypot(DX1, DY1)) # Norma del camp
M[ M == 0] = 1. # Ens evitem que hi hagi divisions per zero
DX1 /= M # Normalitzem
DY1 /= M
#-------------------------------------------------------
# Dibuixem el camp de direccions
# Dibuixem la direccio normalitzada i amb el color de les fletxetes
# la magnitud del camp
p.title("a=%g, b=%g, c= %g" % (a,b,c))
Q = p.quiver(X1, Y1, DX1, DY1, M, pivot='mid', cmap=p.cm.jet)
p.xlabel('x')
p.ylabel('y')
p.legend()
p.grid()
p.xlim(0, xmax)
p.ylim(0, ymax)i obtindrem una gràfica com la següent:
5.3 El model de Lotka-Volterra
Ara anem a fer el model de Lotka-Volterra, és a dir \[x' = u (1-v) y' = a v (v-1)\] on \(a\) és el paràmetre de les equacions. Això ho podem implementar de cop, senzillament canviant la definició del camp:
#bloc6.py
from numpy import *
import pylab as p
from scipy import integrate
# Definim els parametres
a=2
def LotkaVolterra(X, t=0):
""" Retorna les equacions. El t=0 es perque el sistema es autonom """
return array([ X[0]*(1-X[1]),a*X[1]*(X[0]-1.0)])
# Passem a la integracio
t = linspace(0, 40, 1000) # temps integracio
X0 = array([0.1, 0.8]) # condicions inicials: x=0.1, y=0.8
X = integrate.odeint(LotkaVolterra, X0, t)
# Passem a representar-ho:
presa, predador = X.T
#Aixo seria el mateix que si fessim:
#presa=X[:,0]
#predador=X[:,1]
f1 = p.figure() #iniciem una figura
p.plot(t, presa, 'r-', label='presa x (x(0)=%g)' %(presa[0]))
p.plot(t, predador, 'b-', label='predador y (y(0)=%g)' %(predador[0]))
p.grid() # Ens mostra una graella
p.legend(loc='best') # Ens mostra una llegenda al millor lloc
p.xlabel('$t$')
p.ylabel('poblacio')
p.title("a=%g" % (a,))
p.show()Per entendre les oscil·lacions podem mirar al retrat de fase:
#bloc7.py
from numpy import *
import pylab as p
from scipy import integrate
# Definim els parametres
a=0.5
def LotkaVolterra(X, t=0):
""" Retorna les equacions. El t=0 es perque el sistema es autonom """
return array([ X[0]*(1-X[1]),a*X[1]*(X[0]-1.0)])
valors = linspace(1.0, 4.0, 20) # variarem les ci entre 0. i 1.0
vcolors = p.cm.autumn_r(linspace(0.3, 1., len(valors)))
# Per a cada ci triem un color diferent
t = linspace(0, 100, 10000)
ymax = 4 # limits dels eixos
xmax = 4
# Dibuixem el punt fix no trivial
Xinicial=array([1.0,1.0])
p.plot((Xinicial[0],),(Xinicial[1],),'o')
#-------------------------------------------------------
# dibuixem les trajectories que surten dels eixos x
for v, col in zip(valors, vcolors): # Aixo ens fa una doble iteracio
X0 = v*Xinicial # condicio inicial
X = integrate.odeint(LotkaVolterra, X0, t) # no posem l'infodict aqui
p.plot( X[:,0], X[:,1], color=col)
#-------------------------------------------------------
#Definim una xarxa i calculem el camp a la xarxa
num_punts = 30
x = linspace(0, xmax, num_punts)
y = linspace(0, ymax, num_punts)
X1 , Y1 = meshgrid(x, y) # Creem la xarxa
DX1, DY1 = LotkaVolterra([X1, Y1]) # Taxa de creixement a la xarxa
M = (hypot(DX1, DY1)) # Norma del camp
M[ M == 0] = 1. # Ens evitem que hi hagi divisions per zero
DX1 /= M # Normalitzem
DY1 /= M
#-------------------------------------------------------
# Dibuixem el camp de direccions
# Dibuixem la direccio normalitzada i amb el color de les fletxetes
# la magnitud del camp
p.title("a=%g, "% (a,))
Q = p.quiver(X1, Y1, DX1, DY1, M, pivot='mid', cmap=p.cm.jet)
p.xlabel('presa')
p.ylabel('predador')
p.legend()
p.grid()
p.xlim(0, xmax)
p.ylim(0, ymax)
p.show()5.4 Exercici: models més realistes
El model de Lotka-Volterra no és molt realista perquè les oscil·lacions que produeix no són “estructuralment estables” una pertorbació general les destrueix. Un dels models més habituals és el de McArthur-Rozenberg: \[x'=a x \left( 1- \frac{x}{K}\right) - b x y y'= c x y - d y\] Veieu que la següent adimensionalització \[s= a t, \quad \kappa= \frac{c}{d}K, \quad g = \frac{d}{a}\] amb les variables \[u(s)= \frac{c}{d}x(s/a), \quad v(s)= \frac{b}{a} y(s/a)\] produeix el sistema
\[u' = u\left(1-\frac{u}{\kappa}\right)- u v v' = g (u-1) v\] Per aquest model
Trobeu els punts fixos i determineu la seva estabilitat en funció de la diferencial.
Dibuixeu el retrat de fase per algun valor (per exemple, \(\kappa=2\) i \(g=1\)).