#!python
from numpy import *
import pylab as p
# Definition of parameters
a = 1.
b = 0.1
c = 2.0
d = 0.75
def dX_dt(X, t=0):
""" Return the growth rate of fox and rabbit populations. """
return array([ a*X[0] - b*X[0]*X[1] ,
-c*X[1] + d*b*X[0]*X[1] ])
#!python
X_f0 = array([ 0. , 0.])
X_f1 = array([ c/(d*b), a/b])
all(dX_dt(X_f0) == zeros(2) ) and all(dX_dt(X_f1) == zeros(2)) # => True
def d2X_dt2(X, t=0):
""" Return the Jacobian matrix evaluated in X. """
return array([[a -b*X[1], -b*X[0] ],
[b*d*X[1] , -c +b*d*X[0]] ])
A_f0 = d2X_dt2(X_f0) # >>> array([[ 1. , -0. ],
# [ 0. , -1.5]])
#!python
A_f1 = d2X_dt2(X_f1) # >>> array([[ 0. , -2. ],
# [ 0.75, 0. ]])
# whose eigenvalues are +/- sqrt(c*a).j:
lambda1, lambda2 = linalg.eigvals(A_f1) # >>> (1.22474j, -1.22474j)
# They are imaginary numbers. The fox and rabbit populations are periodic as follows from further
# analysis. Their period is given by:
T_f1 = 2*pi/abs(lambda1) # >>> 5.130199
#!python
from scipy import integrate
t = linspace(0, 15, 1000) # time
X0 = array([10, 5]) # initials conditions: 10 rabbits and 5 foxes
X, infodict = integrate.odeint(dX_dt, X0, t, full_output=True)
infodict['message'] # >>> 'Integration successful.'
#!python
rabbits, foxes = X.T
f1 = p.figure()
p.plot(t, rabbits, 'y-', label='Rabbits')
p.plot(t, foxes , 'r-', label='Foxes')
p.plot([], [], ' ', label="Parameters --> a = 1.0,b = 0.1,p = 2.0, q = 0.75")
p.grid()
p.legend(loc='best')
p.xlabel('time')
p.ylabel('population')
p.title('Evolution of fox and rabbit populations')
f1.savefig('rabbits_and_foxes_1.png')
values = linspace(0.3, 0.9, 5) # position of X0 between X_f0 and X_f1
vcolors = p.cm.autumn_r(linspace(0.4, 1., len(values))) # colors for each trajectory
f2 = p.figure()
#-------------------------------------------------------
# plot trajectories
for v, col in zip(values, vcolors):
X0 = v * X_f1 # starting point
X = integrate.odeint( dX_dt, X0, t) # we don't need infodict here
p.plot( X[:,0], X[:,1], lw=3.5*v, color=col, label='X0=(%.f, %.f)' % ( X0[0], X0[1]) )
#-------------------------------------------------------
# define a grid and compute direction at each point
ymax = p.ylim(ymin=0)[1] # get axis limits
xmax = p.xlim(xmin=0)[1]
nb_points = 20
x = linspace(0, xmax, nb_points)
y = linspace(0, ymax, nb_points)
X1 , Y1 = meshgrid(x, y) # create a grid
DX1, DY1 = dX_dt([X1, Y1]) # compute growth rate on the gridt
M = (hypot(DX1, DY1)) # Norm of the growth rate
M[ M == 0] = 1. # Avoid zero division errors
DX1 /= M # Normalize each arrows
DY1 /= M
#-------------------------------------------------------
# Drow direction fields, using matplotlib 's quiver function
# I choose to plot normalized arrows and to use colors to give information on
# the growth speed
p.title('Trajectories and direction fields')
Q = p.quiver(X1, Y1, DX1, DY1, M, pivot='mid', cmap=p.cm.jet)
p.xlabel('Number of rabbits')
p.ylabel('Number of foxes')
p.legend()
p.grid()
p.xlim(0, xmax)
p.ylim(0, ymax)
f2.savefig('rabbits_and_foxes_2.png')
def IF(X):
u, v = X
return u**(c/a) * v * exp( -(b/a)*(d*u+v) )
# We will verify that IF remains constant for different trajectories
for v in values:
X0 = v * X_f1 # starting point
X = integrate.odeint( dX_dt, X0, t)
I = IF(X.T) # compute IF along the trajectory
I_mean = I.mean()
delta = 100 * (I.max()-I.min())/I_mean
print 'X0=(%2.f,%2.f) => I ~ %.1f |delta = %.3G %%' % (X0[0], X0[1], I_mean, delta)
# >>> X0=( 6, 3) => I ~ 20.8 |delta = 6.19E-05 %
# X0=( 9, 4) => I ~ 39.4 |delta = 2.67E-05 %
# X0=(12, 6) => I ~ 55.7 |delta = 1.82E-05 %
# X0=(15, 8) => I ~ 66.8 |delta = 1.12E-05 %
# X0=(18, 9) => I ~ 72.4 |delta = 4.68E-06 %
#!python
from numpy import *
import pylab as p
# Definition of parameters
a = 1.
b = 0.1
c = 1.5
d = 0.75
e= 0.5
f=2.0
g=0.3
def dX_dt(X, t=0):
""" Return the growth rate of fox and rabbit populations. """
return array([ a*X[0] - b*X[0]*X[1] ,
-c*X[1] + d*X[0]*X[1]- e*X[1]*X[2], -f*X[2]+ g*X[1]*X[2] ])
#!python
from scipy import integrate
t = linspace(0, 15, 1000) # time
X0 = array([5, 5, 5]) # initials conditions: 10 rabbits and 5 foxes
X, infodict = integrate.odeint(dX_dt, X0, t, full_output=True)
infodict['message']
#!python
rabbits, foxes, lion = X.T
f1 = p.figure(figsize=(15,5))
p.plot(t, rabbits, 'y-', label='Rabbits')
p.plot(t, foxes , 'r-', label='Foxes')
p.plot(t, lion , 'g-', label='lion')
p.plot([], [], ' ', label="Parameters --> a = 1.0,b = 0.1,p = 2.0, q = 0.75, e=0.5, f=2.0, g=0.1")
p.grid()
p.legend(loc='best')
p.xlabel('time')
p.ylabel('population')
p.title('Evolution of lion, fox and rabbit populations')
f1.savefig('rabbits_and_foxes_lion.png')
values = linspace(0.3, 0.9, 5) # position of X0 between X_f0 and X_f1
vcolors = p.cm.autumn_r(linspace(0.4, 1., len(values))) # colors for each trajectory
X0 = array([5, 5, 5])
f2 = p.figure()
#-------------------------------------------------------
# plot trajectories
for v, col in zip(values, vcolors):
X0 = v * X_f1 # starting point
X = integrate.odeint( dX_dt, X0, t) # we don't need infodict here
p.plot( X[:,0], X[:,1], lw=3.5*v, color=col, label='X0=(%.f, %.f)' % ( X0[0], X0[1]) )
infodict['message']
#-------------------------------------------------------
# define a grid and compute direction at each point
ymax = p.ylim(ymin=0)[1] # get axis limits
xmax = p.xlim(xmin=0)[1]
nb_points = 20
x = linspace(0, xmax, nb_points)
y = linspace(0, ymax, nb_points)
X1 , Y1 = meshgrid(x, y) # create a grid
DX1, DY1 = dX_dt([X1, Y1]) # compute growth rate on the gridt
M = (hypot(DX1, DY1)) # Norm of the growth rate
M[ M == 0] = 1. # Avoid zero division errors
DX1 /= M # Normalize each arrows
DY1 /= M
#-------------------------------------------------------
# Drow direction fields, using matplotlib 's quiver function
# I choose to plot normalized arrows and to use colors to give information on
# the growth speed
p.title('Trajectories and direction fields')
Q = p.quiver(X1, Y1, DX1, DY1, M, pivot='mid', cmap=p.cm.jet)
p.xlabel('Number of rabbits')
p.ylabel('Number of foxes')
p.legend()
p.grid()
p.xlim(0, xmax)
p.ylim(0, ymax)
f2.savefig('rabbits_and_foxes_lion2.png')
# This import registers the 3D projection, but is otherwise unused.
from mpl_toolkits.mplot3d import Axes3D # noqa: F401 unused import
import matplotlib.pyplot as plt
import numpy as np
fig = plt.figure(figsize=(6,6))
ax = fig.gca(projection='3d')
# Make the grid
x, y, z = np.meshgrid(np.arange(-0.8, 1, 0.2),
np.arange(-0.8, 1, 0.2),
np.arange(-0.8, 1, 0.8))
# Make the direction data for the arrows
u = a*X[0] - b*X[0]*X[1]
v = -c*X[1] + d*X[0]*X[1]- e*X[1]*X[2]
w = -f*X[2]+ g*X[1]*X[2]
ax.quiver(x, y, z, u, v, w, length=0.1, normalize=True)
plt.show()
fig.savefig('rabbits_and_foxes_lion2.png')