11# This example describe how to integrate ODEs with scipy.integrate module and how
22# to use the matplotlib module to plot trajectories, directions field and others
33# useful informations.
4- #
5- # == Presentation of the Lokta -Volterra Model ==
6- #
7- #
8- # We will have a look at the Lokta -Volterra model, laso known as the
4+ #
5+ # == Presentation of the Lotka -Volterra Model ==
6+ #
7+ #
8+ # We will have a look at the Lotka -Volterra model, laso known as the
99# predator-prey equations, re a pair of first order, non-linear, differential
1010# equations frequently used to describe the dynamics of biological systems in
1111# which two species interact, one a predator and one its prey. They were proposed
1212# independently by Alfred J. Lotka in 1925 and Vito Volterra in 1926 :
1313# du/dt = a*u - b*u*v
14- # dv/dt = -c*v + d*b*u*v
15- #
14+ # dv/dt = -c*v + d*b*u*v
15+ #
1616# with the following notations :
17- #
17+ #
1818# * u : number of prey (for example rabbits)
19- #
20- # * v : number of predators (for example foxes)
21- #
22- # * a, b, c, d are constant parameters defining the behavior of the population :
23- #
19+ #
20+ # * v : number of predators (for example foxes)
21+ #
22+ # * a, b, c, d are constant parameters defining the behavior of the population :
23+ #
2424# + a is the natural growing rate of rabbits, when there's no foxes
25- #
25+ #
2626# + b is the natural dying rate of rabbits, due to predation
27- #
27+ #
2828# + c is the natural dying rate of foxes, when there's no rabbits
29- #
29+ #
3030# + d is the factor descibing how many catched rabbits let create new rabbits
31- #
32- #
31+ #
32+ #
3333# We will use X=[u, v] to describe the state of both populations.
3434# Definition of the equations:
35- #
35+ #
3636from numpy import *
3737import pylab as p
3838
39- # Parameters
39+ # Parameters
4040a = 1.
4141b = 0.1
4242c = 1.5
4343d = 0.75
4444
4545def dX_dt (X , t = 0 ):
4646 """ Return the growth rate of foxes and rabbits populations. """
47- return array ([ a * X [0 ] - b * X [0 ]* X [1 ] ,
47+ return array ([ a * X [0 ] - b * X [0 ]* X [1 ] ,
4848 - c * X [1 ] + d * b * X [0 ]* X [1 ] ])
49- #
49+ #
5050# === Population equilibrium ===
51- #
52- # Before using scipy to integrate this system, we will have a closer look on
51+ #
52+ # Before using scipy to integrate this system, we will have a closer look on
5353# position equilibrium. Equilibrium occurs when the growth rate is equal to 0.
5454# This gives two fixed points:
55- #
55+ #
5656X_f0 = array ([0. , 0. ])
5757X_f1 = array ([ c / (d * b ), a / b ])
58- all (dX_dt (X_f0 ) == zeros (2 ) ) and all (dX_dt (X_f1 ) == zeros (2 )) # => True
59- #
58+ all (dX_dt (X_f0 ) == zeros (2 ) ) and all (dX_dt (X_f1 ) == zeros (2 )) # => True
59+ #
6060# === Stability of the fixed points ===
6161# Near theses two points, the system can be linearized :
6262# dX_dt = A_f*X where A is the Jacobian matrix evaluated at the corresponding point.
6363# We have to define the Jacobian matrix:
64- #
64+ #
6565def d2X_dt2 (X , t = 0 ):
6666 """ Return the jacobian matrix evaluated in X. """
6767 return array ([[a - b * X [1 ], - b * X [0 ] ],
68- [b * d * X [1 ] , - c + b * d * X [0 ]] ])
69- #
68+ [b * d * X [1 ] , - c + b * d * X [0 ]] ])
69+ #
7070# So, near X_f0, which represents the extinction of both species, we have:
7171# A_f0 = d2X_dt2(X_f0) # >>> array([[ 1. , -0. ],
7272# # [ 0. , -1.5]])
73- #
73+ #
7474# Near X_f0, the number of rabbits increase and the population of foxes decrease.
7575# X_f0 is a [http://en.wikipedia.org/wiki/Saddle_point saddle point].
76- #
76+ #
7777# Near X_f1, we have :
7878A_f1 = d2X_dt2 (X_f1 ) # >>> array([[ 0. , -2. ],
7979 # [ 0.75, 0. ]])
@@ -84,25 +84,25 @@ def d2X_dt2(X, t=0):
8484# They are imaginary number, so the fox and rabbit populations are periodic and
8585# their period is given by :
8686T_f1 = 2 * pi / abs (lambda1 ) # >>> 5.130199
87- #
87+ #
8888# == Integrating the ODE using scipy.integate ==
89- #
89+ #
9090# Now we will use the scipy.integrate module to integrate the ODE.
9191# This module offers a method named odeint, very easy to use to integrade ODE:
92- #
92+ #
9393from scipy import integrate
9494
9595t = linspace (0 , 15 , 1000 ) # time
96- X0 = array ([10 , 5 ]) # initials conditions: 10 rabbits and 5 foxes
96+ X0 = array ([10 , 5 ]) # initials conditions: 10 rabbits and 5 foxes
9797
9898X , infodict = integrate .odeint (dX_dt , X0 , t , full_output = True )
9999infodict ['message' ] # >>> 'Integration successful.'
100- #
100+ #
101101# `infodict` is optionnal, and you can omit the `full_output` argument if you don't want it.
102102# type "info(odeint)" if you want more information about odeint inputs and outputs.
103- #
103+ #
104104# We will use matplotlib to plot the evolution of both populations:
105- #
105+ #
106106rabbits , foxes = X .T
107107
108108f1 = p .figure ()
@@ -114,46 +114,46 @@ def d2X_dt2(X, t=0):
114114p .ylabel ('population' )
115115p .title ('Evolution of fox and rabbit populations' )
116116f1 .savefig ('rabbits_and_foxes_1.png' )
117- #
118- #
117+ #
118+ #
119119# The populations are indeed periodic, and their period is near to the T_f1 we calculated.
120- #
120+ #
121121# == Plotting directions field and trajectories in the phase plane ==
122- #
122+ #
123123# We will plot some trajectories in a phase plane for different starting
124124# points between X__f0 and X_f1.
125- #
125+ #
126126# We will ue matplotlib's colormap to define colors for the trajectories.
127127# These colormaps are very useful to make nice plots.
128128# Have a look on ShowColormaps if you want more information.
129- #
129+ #
130130values = linspace (0.3 , 0.9 , 5 ) # position of X0 between X_f0 and X_f1
131131vcolors = p .cm .Greens (linspace (0.3 , 1. , len (values ))) # colors for each trajectory
132132
133133f2 = p .figure ()
134134
135135#-------------------------------------------------------
136136# plot trajectories
137- for v , col in zip (values , vcolors ):
137+ for v , col in zip (values , vcolors ):
138138 X0 = v * X_f1 # starting point
139139 X = integrate .odeint ( dX_dt , X0 , t ) # we don't need infodict here
140140 p .plot ( X [:,0 ], X [:,1 ], lw = 3.5 * v , color = col , label = 'X0=(%.f, %.f)' % ( X0 [0 ], X0 [1 ]) )
141141
142142#-------------------------------------------------------
143143# define a grid and compute direction at each point
144144ymax = p .ylim (ymin = 0 )[1 ] # get axis limits
145- xmax = p .xlim (xmin = 0 )[1 ]
146- nb_points = 20
145+ xmax = p .xlim (xmin = 0 )[1 ]
146+ nb_points = 20
147147
148148x = linspace (0 , xmax , nb_points )
149149y = linspace (0 , ymax , nb_points )
150150
151151X1 , Y1 = meshgrid (x , y ) # create a grid
152152DX1 , DY1 = dX_dt ([X1 , Y1 ]) # compute growth rate on the gridt
153- M = (hypot (DX1 , DY1 )) # Norm of the growth rate
154- M [ M == 0 ] = 1. # Avoid zero division errors
153+ M = (hypot (DX1 , DY1 )) # Norm of the growth rate
154+ M [ M == 0 ] = 1. # Avoid zero division errors
155155DX1 /= M # Normalize each arrows
156- DY1 /= M
156+ DY1 /= M
157157
158158#-------------------------------------------------------
159159# Drow direction fields, using matplotlib 's quiver function
@@ -168,18 +168,18 @@ def d2X_dt2(X, t=0):
168168p .xlim (0 , xmax )
169169p .ylim (0 , ymax )
170170f2 .savefig ('rabbits_and_foxes_2.png' )
171- #
172- #
171+ #
172+ #
173173# We can see on this graph that an intervention on fox or rabbit populations can
174174# have non intuitive effects. If, in order to decrease the number of rabbits,
175175# we introduce foxes, this can lead to an increase of rabbits in the long run,
176176# if that intervention happens at a bad moment.
177- #
178- #
177+ #
178+ #
179179# == Plotting contours ==
180- #
180+ #
181181# We can verify that the fonction IF defined below remains constant along a trajectory:
182- #
182+ #
183183def IF (X ):
184184 u , v = X
185185 return u ** (c / a ) * v * exp ( - (b / a )* (d * u + v ) )
@@ -189,9 +189,9 @@ def IF2(X):
189189 return u ** (c / a ) * v * exp ( - (b / a )* (d * u + v ) )
190190
191191# We will verify that IF remains constant for differents trajectories
192- for v in values :
192+ for v in values :
193193 X0 = v * X_f1 # starting point
194- X = integrate .odeint ( dX_dt , X0 , t )
194+ X = integrate .odeint ( dX_dt , X0 , t )
195195 I = IF (X .T ) # compute IF along the trajectory
196196 I_mean = I .mean ()
197197 delta = 100 * (I .max ()- I .min ())/ I_mean
@@ -202,15 +202,15 @@ def IF2(X):
202202# X0=(12, 6) => I ~ 55.7 |delta = 1.82E-05 %
203203# X0=(15, 8) => I ~ 66.8 |delta = 1.12E-05 %
204204# X0=(18, 9) => I ~ 72.4 |delta = 4.68E-06 %
205- #
205+ #
206206# Potting iso-contours of IF can be a good representation of trajectories,
207207# without having to integrate the ODE
208- #
208+ #
209209#-------------------------------------------------------
210210# plot iso contours
211- nb_points = 80 # grid size
211+ nb_points = 80 # grid size
212212
213- x = linspace (0 , xmax , nb_points )
213+ x = linspace (0 , xmax , nb_points )
214214y = linspace (0 , ymax , nb_points )
215215
216216X2 , Y2 = meshgrid (x , y ) # create the grid
@@ -228,6 +228,6 @@ def IF2(X):
228228p .title ('IF contours' )
229229f3 .savefig ('rabbits_and_foxes_3.png' )
230230p .show ()
231- #
232- #
231+ #
232+ #
233233# # vim: set et sts=4 sw=4:
0 commit comments