1
1
# This example describe how to integrate ODEs with scipy.integrate module and how
2
2
# to use the matplotlib module to plot trajectories, directions field and others
3
3
# 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
9
9
# predator-prey equations, re a pair of first order, non-linear, differential
10
10
# equations frequently used to describe the dynamics of biological systems in
11
11
# which two species interact, one a predator and one its prey. They were proposed
12
12
# independently by Alfred J. Lotka in 1925 and Vito Volterra in 1926 :
13
13
# 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
+ #
16
16
# with the following notations :
17
- #
17
+ #
18
18
# * 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
+ #
24
24
# + a is the natural growing rate of rabbits, when there's no foxes
25
- #
25
+ #
26
26
# + b is the natural dying rate of rabbits, due to predation
27
- #
27
+ #
28
28
# + c is the natural dying rate of foxes, when there's no rabbits
29
- #
29
+ #
30
30
# + d is the factor descibing how many catched rabbits let create new rabbits
31
- #
32
- #
31
+ #
32
+ #
33
33
# We will use X=[u, v] to describe the state of both populations.
34
34
# Definition of the equations:
35
- #
35
+ #
36
36
from numpy import *
37
37
import pylab as p
38
38
39
- # Parameters
39
+ # Parameters
40
40
a = 1.
41
41
b = 0.1
42
42
c = 1.5
43
43
d = 0.75
44
44
45
45
def dX_dt (X , t = 0 ):
46
46
""" 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 ] ,
48
48
- c * X [1 ] + d * b * X [0 ]* X [1 ] ])
49
- #
49
+ #
50
50
# === 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
53
53
# position equilibrium. Equilibrium occurs when the growth rate is equal to 0.
54
54
# This gives two fixed points:
55
- #
55
+ #
56
56
X_f0 = array ([0. , 0. ])
57
57
X_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
+ #
60
60
# === Stability of the fixed points ===
61
61
# Near theses two points, the system can be linearized :
62
62
# dX_dt = A_f*X where A is the Jacobian matrix evaluated at the corresponding point.
63
63
# We have to define the Jacobian matrix:
64
- #
64
+ #
65
65
def d2X_dt2 (X , t = 0 ):
66
66
""" Return the jacobian matrix evaluated in X. """
67
67
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
+ #
70
70
# So, near X_f0, which represents the extinction of both species, we have:
71
71
# A_f0 = d2X_dt2(X_f0) # >>> array([[ 1. , -0. ],
72
72
# # [ 0. , -1.5]])
73
- #
73
+ #
74
74
# Near X_f0, the number of rabbits increase and the population of foxes decrease.
75
75
# X_f0 is a [http://en.wikipedia.org/wiki/Saddle_point saddle point].
76
- #
76
+ #
77
77
# Near X_f1, we have :
78
78
A_f1 = d2X_dt2 (X_f1 ) # >>> array([[ 0. , -2. ],
79
79
# [ 0.75, 0. ]])
@@ -84,25 +84,25 @@ def d2X_dt2(X, t=0):
84
84
# They are imaginary number, so the fox and rabbit populations are periodic and
85
85
# their period is given by :
86
86
T_f1 = 2 * pi / abs (lambda1 ) # >>> 5.130199
87
- #
87
+ #
88
88
# == Integrating the ODE using scipy.integate ==
89
- #
89
+ #
90
90
# Now we will use the scipy.integrate module to integrate the ODE.
91
91
# This module offers a method named odeint, very easy to use to integrade ODE:
92
- #
92
+ #
93
93
from scipy import integrate
94
94
95
95
t = 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
97
97
98
98
X , infodict = integrate .odeint (dX_dt , X0 , t , full_output = True )
99
99
infodict ['message' ] # >>> 'Integration successful.'
100
- #
100
+ #
101
101
# `infodict` is optionnal, and you can omit the `full_output` argument if you don't want it.
102
102
# type "info(odeint)" if you want more information about odeint inputs and outputs.
103
- #
103
+ #
104
104
# We will use matplotlib to plot the evolution of both populations:
105
- #
105
+ #
106
106
rabbits , foxes = X .T
107
107
108
108
f1 = p .figure ()
@@ -114,46 +114,46 @@ def d2X_dt2(X, t=0):
114
114
p .ylabel ('population' )
115
115
p .title ('Evolution of fox and rabbit populations' )
116
116
f1 .savefig ('rabbits_and_foxes_1.png' )
117
- #
118
- #
117
+ #
118
+ #
119
119
# The populations are indeed periodic, and their period is near to the T_f1 we calculated.
120
- #
120
+ #
121
121
# == Plotting directions field and trajectories in the phase plane ==
122
- #
122
+ #
123
123
# We will plot some trajectories in a phase plane for different starting
124
124
# points between X__f0 and X_f1.
125
- #
125
+ #
126
126
# We will ue matplotlib's colormap to define colors for the trajectories.
127
127
# These colormaps are very useful to make nice plots.
128
128
# Have a look on ShowColormaps if you want more information.
129
- #
129
+ #
130
130
values = linspace (0.3 , 0.9 , 5 ) # position of X0 between X_f0 and X_f1
131
131
vcolors = p .cm .Greens (linspace (0.3 , 1. , len (values ))) # colors for each trajectory
132
132
133
133
f2 = p .figure ()
134
134
135
135
#-------------------------------------------------------
136
136
# plot trajectories
137
- for v , col in zip (values , vcolors ):
137
+ for v , col in zip (values , vcolors ):
138
138
X0 = v * X_f1 # starting point
139
139
X = integrate .odeint ( dX_dt , X0 , t ) # we don't need infodict here
140
140
p .plot ( X [:,0 ], X [:,1 ], lw = 3.5 * v , color = col , label = 'X0=(%.f, %.f)' % ( X0 [0 ], X0 [1 ]) )
141
141
142
142
#-------------------------------------------------------
143
143
# define a grid and compute direction at each point
144
144
ymax = 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
147
147
148
148
x = linspace (0 , xmax , nb_points )
149
149
y = linspace (0 , ymax , nb_points )
150
150
151
151
X1 , Y1 = meshgrid (x , y ) # create a grid
152
152
DX1 , 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
155
155
DX1 /= M # Normalize each arrows
156
- DY1 /= M
156
+ DY1 /= M
157
157
158
158
#-------------------------------------------------------
159
159
# Drow direction fields, using matplotlib 's quiver function
@@ -168,18 +168,18 @@ def d2X_dt2(X, t=0):
168
168
p .xlim (0 , xmax )
169
169
p .ylim (0 , ymax )
170
170
f2 .savefig ('rabbits_and_foxes_2.png' )
171
- #
172
- #
171
+ #
172
+ #
173
173
# We can see on this graph that an intervention on fox or rabbit populations can
174
174
# have non intuitive effects. If, in order to decrease the number of rabbits,
175
175
# we introduce foxes, this can lead to an increase of rabbits in the long run,
176
176
# if that intervention happens at a bad moment.
177
- #
178
- #
177
+ #
178
+ #
179
179
# == Plotting contours ==
180
- #
180
+ #
181
181
# We can verify that the fonction IF defined below remains constant along a trajectory:
182
- #
182
+ #
183
183
def IF (X ):
184
184
u , v = X
185
185
return u ** (c / a ) * v * exp ( - (b / a )* (d * u + v ) )
@@ -189,9 +189,9 @@ def IF2(X):
189
189
return u ** (c / a ) * v * exp ( - (b / a )* (d * u + v ) )
190
190
191
191
# We will verify that IF remains constant for differents trajectories
192
- for v in values :
192
+ for v in values :
193
193
X0 = v * X_f1 # starting point
194
- X = integrate .odeint ( dX_dt , X0 , t )
194
+ X = integrate .odeint ( dX_dt , X0 , t )
195
195
I = IF (X .T ) # compute IF along the trajectory
196
196
I_mean = I .mean ()
197
197
delta = 100 * (I .max ()- I .min ())/ I_mean
@@ -202,15 +202,15 @@ def IF2(X):
202
202
# X0=(12, 6) => I ~ 55.7 |delta = 1.82E-05 %
203
203
# X0=(15, 8) => I ~ 66.8 |delta = 1.12E-05 %
204
204
# X0=(18, 9) => I ~ 72.4 |delta = 4.68E-06 %
205
- #
205
+ #
206
206
# Potting iso-contours of IF can be a good representation of trajectories,
207
207
# without having to integrate the ODE
208
- #
208
+ #
209
209
#-------------------------------------------------------
210
210
# plot iso contours
211
- nb_points = 80 # grid size
211
+ nb_points = 80 # grid size
212
212
213
- x = linspace (0 , xmax , nb_points )
213
+ x = linspace (0 , xmax , nb_points )
214
214
y = linspace (0 , ymax , nb_points )
215
215
216
216
X2 , Y2 = meshgrid (x , y ) # create the grid
@@ -228,6 +228,6 @@ def IF2(X):
228
228
p .title ('IF contours' )
229
229
f3 .savefig ('rabbits_and_foxes_3.png' )
230
230
p .show ()
231
- #
232
- #
231
+ #
232
+ #
233
233
# # vim: set et sts=4 sw=4:
0 commit comments