Commit e9314d2a authored by Romain Casati's avatar Romain Casati
Browse files

Adding Lotka-Volterra example.

parent 697623a2
Pipeline #413773 passed with stages
in 3 minutes and 47 seconds
{"cells":[{"metadata":{},"cell_type":"markdown","source":"Taken from the [scipy cookbook](https://scipy-cookbook.readthedocs.io/items/LoktaVolterraTutorial.html)"},{"metadata":{},"cell_type":"markdown","source":"# Matplotlib: lotka volterra tutorial\n\n\nThis example describes how to integrate ODEs with the scipy.integrate\nmodule, and how to use the matplotlib module to plot trajectories,\ndirection fields and other information.\n\nYou can get the source code for this tutorial here:\n[tutorial_lokta-voltera_v4.py](files/attachments/LoktaVolterraTutorial/tutorial_lokta-voltera_v4.py).\n\nPresentation of the Lotka-Volterra Model\n----------------------------------------\n\nWe will have a look at the Lotka-Volterra model, also known as the\npredator-prey equations, which is a pair of first order, non-linear,\ndifferential equations frequently used to describe the dynamics of\nbiological systems in which two species interact, one a predator and the\nother its prey. The model was proposed independently by Alfred J. Lotka\nin 1925 and Vito Volterra in 1926, and can be described by"},{"metadata":{},"cell_type":"markdown","source":"$$ du/dt = a*u - b*u*v$$\n$$dv/dt = -c*v + d*b*u*v$$"},{"metadata":{},"cell_type":"markdown","source":"with the following notations:\n\n* u: number of preys (for example, rabbits)\n\n* v: number of predators (for example, foxes) \n \n* a, b, c, d are constant parameters defining the behavior of the population: \n\n + a is the natural growing rate of rabbits, when there's no fox\n\n + b is the natural dying rate of rabbits, due to predation\n\n + c is the natural dying rate of fox, when there's no rabbit\n\n + d is the factor describing how many caught rabbits let create a new fox"},{"metadata":{},"cell_type":"markdown","source":"We will use X=[u, v] to describe the state of both populations.\n\nDefinition of the equations:"},{"metadata":{"trusted":true},"cell_type":"code","source":"from numpy import *\nimport matplotlib.pyplot as p\n# Definition of parameters\na = 1.\nb = 0.1\nc = 1.5\nd = 0.75\ndef dX_dt(X, t=0):\n \"\"\" Return the growth rate of fox and rabbit populations. \"\"\"\n return array([ a*X[0] - b*X[0]*X[1] ,\n -c*X[1] + d*b*X[0]*X[1] ])","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"### Population equilibrium\n\nBefore using !SciPy to integrate this system, we will have a closer look\nat position equilibrium. Equilibrium occurs when the growth rate is\nequal to 0. This gives two fixed points:"},{"metadata":{"trusted":true},"cell_type":"code","source":"X_f0 = array([ 0. , 0.])\nX_f1 = array([ c/(d*b), a/b])\nall(dX_dt(X_f0) == zeros(2) ) and all(dX_dt(X_f1) == zeros(2)) # => True","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"### Stability of the fixed points\n\nNear these two points, the system can be linearized: dX\\_dt = A\\_f\\*X\nwhere A is the Jacobian matrix evaluated at the corresponding point. We\nhave to define the Jacobian matrix:"},{"metadata":{"trusted":true},"cell_type":"code","source":"def d2X_dt2(X, t=0):\n \"\"\" Return the Jacobian matrix evaluated in X. \"\"\"\n return array([[a -b*X[1], -b*X[0] ],\n [b*d*X[1] , -c +b*d*X[0]] ])","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"So near X\\_f0, which represents the extinction of both species, we have:"},{"metadata":{"trusted":true},"cell_type":"code","source":"A_f0 = d2X_dt2(X_f0) # >>> array([[ 1. , -0. ],\n # [ 0. , -1.5]])","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"Near X\\_f0, the number of rabbits increase and the population of foxes\ndecrease. The origin is therefore a [saddle\npoint](http://en.wikipedia.org/wiki/Saddle_point).\n\nNear X\\_f1, we have:"},{"metadata":{"trusted":true},"cell_type":"code","source":"A_f1 = d2X_dt2(X_f1) # >>> array([[ 0. , -2. ],\n # [ 0.75, 0. ]])\n# whose eigenvalues are +/- sqrt(c*a).j:\nlambda1, lambda2 = linalg.eigvals(A_f1) # >>> (1.22474j, -1.22474j)\n# They are imaginary numbers. The fox and rabbit populations are periodic as follows from further\n# analysis. Their period is given by:\nT_f1 = 2*pi/abs(lambda1) # >>> 5.130199","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"Integrating the ODE using scipy.integrate\n-----------------------------------------\n\nNow we will use the scipy.integrate module to integrate the ODEs. This\nmodule offers a method named odeint, which is very easy to use to\nintegrate ODEs:"},{"metadata":{"trusted":true},"cell_type":"code","source":"from scipy import integrate\nt = linspace(0, 15, 1000) # time\nX0 = array([10, 5]) # initials conditions: 10 rabbits and 5 foxes\nX = integrate.odeint(dX_dt, X0, t)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"We can now use Matplotlib to plot the evolution of both populations:"},{"metadata":{"trusted":true},"cell_type":"code","source":"rabbits, foxes = X.T\nf1 = p.figure()\np.plot(t, rabbits, 'r-', label='Rabbits')\np.plot(t, foxes , 'b-', label='Foxes')\np.grid()\np.legend(loc='best')\np.xlabel('time')\np.ylabel('population')\np.title('Evolution of fox and rabbit populations')\np.show()","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"![](files/attachments/LoktaVolterraTutorial/rabbits_and_foxes_1v2.png)\n\nThe populations are indeed periodic, and their period is close to the\nvalue T\\_f1 that we computed.\n\nPlotting direction fields and trajectories in the phase plane\n-------------------------------------------------------------\n\nWe will plot some trajectories in a phase plane for different starting\npoints between X\\_f0 and X\\_f1.\n\nWe will use Matplotlib's colormap to define colors for the trajectories.\nThese colormaps are very useful to make nice plots. Have a look at\n[ShowColormaps](http://www.scipy.org/Cookbook/Matplotlib/Show_colormaps)\nif you want more information."},{"metadata":{"trusted":true},"cell_type":"code","source":"values = linspace(0.3, 0.9, 5) # position of X0 between X_f0 and X_f1\nvcolors = p.cm.autumn_r(linspace(0.3, 1., len(values))) # colors for each trajectory\n\nf2 = p.figure()\n\n#-------------------------------------------------------\n# plot trajectories\nfor v, col in zip(values, vcolors): \n X0 = v * X_f1 # starting point\n X = integrate.odeint( dX_dt, X0, t) # we don't need infodict here\n p.plot( X[:,0], X[:,1], lw=3.5*v, color=col, label='X0=(%.f, %.f)' % ( X0[0], X0[1]) )\n\n#-------------------------------------------------------\n# define a grid and compute direction at each point\nymax = p.ylim(ymin=0)[1] # get axis limits\nxmax = p.xlim(xmin=0)[1] \nnb_points = 20 \n\nx = linspace(0, xmax, nb_points)\ny = linspace(0, ymax, nb_points)\n\nX1 , Y1 = meshgrid(x, y) # create a grid\nDX1, DY1 = dX_dt([X1, Y1]) # compute growth rate on the gridt\nM = (hypot(DX1, DY1)) # Norm of the growth rate \nM[ M == 0] = 1. # Avoid zero division errors \nDX1 /= M # Normalize each arrows\nDY1 /= M \n\n#-------------------------------------------------------\n# Drow direction fields, using matplotlib 's quiver function\n# I choose to plot normalized arrows and to use colors to give information on\n# the growth speed\np.title('Trajectories and direction fields')\nQ = p.quiver(X1, Y1, DX1, DY1, M, pivot='mid', cmap=p.cm.jet)\np.xlabel('Number of rabbits')\np.ylabel('Number of foxes')\np.legend()\np.grid()\np.xlim(0, xmax)\np.ylim(0, ymax)\np.show()","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"![](files/attachments/LoktaVolterraTutorial/rabbits_and_foxes_2v3.png)\n\nThis graph shows us that changing either the fox or the rabbit\npopulation can have an unintuitive effect. If, in order to decrease the\nnumber of rabbits, we introduce foxes, this can lead to an increase of\nrabbits in the long run, depending on the time of intervention.\n\nPlotting contours\n-----------------\n\nWe can verify that the function IF defined below remains constant along\na trajectory:"},{"metadata":{"trusted":true},"cell_type":"code","source":"def IF(X):\n u, v = X\n return u**(c/a) * v * exp( -(b/a)*(d*u+v) )\n# We will verify that IF remains constant for different trajectories\nfor v in values:\n X0 = v * X_f1 # starting point\n X = integrate.odeint( dX_dt, X0, t)\n I = IF(X.T) # compute IF along the trajectory\n I_mean = I.mean()\n delta = 100 * (I.max()-I.min())/I_mean\n print('X0=(%2.f,%2.f) => I ~ %.1f |delta = %.3G %%' % (X0[0], X0[1], I_mean, delta))\n# >>> X0=( 6, 3) => I ~ 20.8 |delta = 6.19E-05 %\n# X0=( 9, 4) => I ~ 39.4 |delta = 2.67E-05 %\n# X0=(12, 6) => I ~ 55.7 |delta = 1.82E-05 %\n# X0=(15, 8) => I ~ 66.8 |delta = 1.12E-05 %\n# X0=(18, 9) => I ~ 72.4 |delta = 4.68E-06 %","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"Plotting iso-contours of IF can be a good representation of\ntrajectories, without having to integrate the ODE"},{"metadata":{"trusted":true},"cell_type":"code","source":"# plot iso contours\nnb_points = 80 # grid size\nx = linspace(0, xmax, nb_points)\ny = linspace(0, ymax, nb_points)\nX2 , Y2 = meshgrid(x, y) # create the grid\nZ2 = IF([X2, Y2]) # compute IF on each point\nf3 = p.figure()\nCS = p.contourf(X2, Y2, Z2, cmap=p.cm.Purples_r, alpha=0.5)\nCS2 = p.contour(X2, Y2, Z2, colors='black', linewidths=2. )\np.clabel(CS2, inline=1, fontsize=16, fmt='%.f')\np.grid()\np.xlabel('Number of rabbits')\np.ylabel('Number of foxes')\np.ylim(1, ymax)\np.xlim(1, xmax)\np.title('IF contours')\nf3.savefig('rabbits_and_foxes_3.png')\np.show()","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"![](files/attachments/LoktaVolterraTutorial/rabbits_and_foxes_3v2.png)"}],"metadata":{"kernelspec":{"display_name":"Python 2","language":"python","name":"python2"},"language_info":{"codemirror_mode":{"name":"ipython","version":2},"file_extension":".py","mimetype":"text/x-python","name":"python","nbconvert_exporter":"python","pygments_lexer":"ipython2","version":"2.7.13"}},"nbformat":4,"nbformat_minor":2}
\ No newline at end of file
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment