Modèle Proies-Prédateurs de Lotka-Volterra

Université de Nice, Département de Maths

Le modèle que nous étudions a été proposé par Volterra (et indépendamment par Lotka) en 1926 dans un ouvrage intitulé ”Théorie mathématique de la lutte pour la vie” qui est probablement le premier traité d’écologie mathématique. 

Volterra avait été consulté par le responsable de la pêche italienne `a Trieste qui avait remarqué que, juste après la première guerre mondiale (période durant laquelle la pêche avait été nettement réduite) la proportion de requins et autres prédateurs impropres `a la consommation que l’on pêchait parmi les poissons consommables était nettement supérieure `a ce qu’elle était avant guerre et `a ce qu’elle redevint ensuite.

https://math.unice.fr/~diener/MAB06/LotVolt.pdf

# This example describe how to integrate ODEs with scipy.integrate module, and how

# to use the matplotlib module to plot trajectories, direction fields and other

# useful information.

# 

# == Presentation of the Lokta-Volterra Model ==

# 

# We will have a look at the Lokta-Volterra model, also known as the

# predator-prey equations, which are a pair of first order, non-linear, differential

# equations frequently used to describe the dynamics of biological systems in

# which two species interact, one a predator and one its prey. They were proposed

# independently by Alfred J. Lotka in 1925 and Vito Volterra in 1926:

# du/dt =  a*u -   b*u*v

# dv/dt = -c*v + d*b*u*v 

# 

# with the following notations:

# 

# *  u: number of preys (for example, rabbits)

# 

# *  v: number of predators (for example, foxes)  

#   

# * a, b, c, d are constant parameters defining the behavior of the population:    

# 

#   + a is the natural growing rate of rabbits, when there's no fox

# 

#   + b is the natural dying rate of rabbits, due to predation

# 

#   + c is the natural dying rate of fox, when there's no rabbit

# 

#   + d is the factor describing how many caught rabbits let create a new fox

# 

# We will use X=[u, v] to describe the state of both populations.

# 

# Definition of the equations:

# 

from numpy import *

import pylab as p

# Definition of parameters 

a = 1.

b = 0.1

c = 1.5

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] ])

# 

# === Population equilibrium ===

# 

# Before using !SciPy to integrate this system, we will have a closer look on 

# position equilibrium. Equilibrium occurs when the growth rate is equal to 0.

# This gives two fixed points:

# 

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 

# 

# === Stability of the fixed points ===

# Near theses two points, the system can be linearized:

# dX_dt = A_f*X where A is the Jacobian matrix evaluated at the corresponding point.

# We have to define the Jacobian matrix:

# 

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]] ])  

# 

# So, near X_f0, which represents the extinction of both species, we have:

# A_f0 = d2X_dt2(X_f0)                    # >>> array([[ 1. , -0. ],

#                                         #            [ 0. , -1.5]])

# 

# Near X_f0, the number of rabbits increase and the population of foxes decrease.

# The origin is a [http://en.wikipedia.org/wiki/Saddle_point saddle point].

# 

# Near X_f1, we have:

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 number, so the fox and rabbit populations are periodic and

# their period is given by:

T_f1 = 2*pi/abs(lambda1)                # >>> 5.130199

#         

# == Integrating the ODE using scipy.integate ==

# 

# Now we will use the scipy.integrate module to integrate the ODEs.

# This module offers a method named odeint, very easy to use to integrate ODEs:

# 

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.'

# 

# `infodict` is optional, and you can omit the `full_output` argument if you don't want it.

# Type "info(odeint)" if you want more information about odeint inputs and outputs.

# 

# We can now use Matplotlib to plot the evolution of both populations:

# 

rabbits, foxes = X.T

f1 = p.figure()

p.plot(t, rabbits, 'r-', label='Rabbits')

p.plot(t, foxes  , 'b-', label='Foxes')

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')

# 

# 

# The populations are indeed periodic, and their period is near to the T_f1 we calculated.

# 

# == Plotting direction fields and trajectories in the phase plane ==

# 

# We will plot some trajectories in a phase plane for different starting

# points between X__f0 and X_f1.

# 

# We will use matplotlib's colormap to define colors for the trajectories.

# These colormaps are very useful to make nice plots.

# Have a look at [http://www.scipy.org/Cookbook/Matplotlib/Show_colormaps ShowColormaps] if you want more information.

# 

values  = linspace(0.3, 0.9, 5)                          # position of X0 between X_f0 and X_f1

vcolors = p.cm.autumn_r(linspace(0.3, 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')

# 

# 

# We can see on this graph that an intervention on fox or rabbit populations can

# have non intuitive effects. If, in order to decrease the number of rabbits,

# we introduce foxes, this can lead to an increase of rabbits in the long run,

# if that intervention happens at a bad moment.

# 

# 

# == Plotting contours ==

# 

# We can verify that the function IF defined below remains constant along a trajectory:

# 

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 %

# 

# Potting iso-contours of IF can be a good representation of trajectories,

# without having to integrate the ODE

# 

#-------------------------------------------------------

# plot iso contours

nb_points = 80                              # grid size 

x = linspace(0, xmax, nb_points)    

y = linspace(0, ymax, nb_points)

X2 , Y2  = meshgrid(x, y)                   # create the grid

Z2 = IF([X2, Y2])                           # compute IF on each point

f3 = p.figure()

CS = p.contourf(X2, Y2, Z2, cmap=p.cm.Purples_r, alpha=0.5)

CS2 = p.contour(X2, Y2, Z2, colors='black', linewidths=2. )

p.clabel(CS2, inline=1, fontsize=16, fmt='%.f')

p.grid()

p.xlabel('Number of rabbits')

p.ylabel('Number of foxes')

p.ylim(1, ymax)

p.xlim(1, xmax)

p.title('IF contours')

f3.savefig('rabbits_and_foxes_3.png')

p.show()

# 

# 

# # vim: set et sts=4 sw=4: