Plotting trajectories and direction fields for a system of ODEs in Python

The aim of this post is to guide the reader through plotting trajectories and direction fields for a system of equations in Python. This is useful when investigating the equilibria and stability of the system, and to facilitate in understanding the general behavior of a system under study. I will use a system of predator-prey equations, that my most devoted online readers are already familiar with from my previous posts on identifying equilibria and stability, and on nondimensionalization. Specifically, I’ll be using the Lotka-Volterra set of equations with Holling’s Type II functional response:

\frac{\mathrm{d} x}{\mathrm{d} t}=bx\left ( 1-\frac{x}{K} \right )-\frac{axy}{1+ahx}

\frac{\mathrm{d} y}{\mathrm{d} t}=\frac{caxy}{1+ahx}-dy

where:

x: prey abundance

y: predator abundance

b: prey growth rate

d: predator death rate

c: rate with which consumed prey is converted to predator

a: rate with which prey is killed by a predator per unit of time

K: prey carrying capacity given the prey’s environmental conditions

h: handling time

This system has 3 equilibria: when both species are dead (0,0), when predators are dead and the prey grows to its carrying capacity (K,0) and a non-trivial equilibrium where both species coexist and is generally more interesting, given by:

y^*=\frac{b}{a}(1+ahx^*)\left(1-\frac{x^*}{K} \right)

x^*=\frac{d}{a(c-dh)}

The following code should produce both trajectories and direction fields for this system of ODEs (python virtuosos please excuse the extensive commenting, I try to comment as much as possible for people new to python):

import numpy as np
from matplotlib import pyplot as plt
from scipy import integrate

# I'm using this style for a pretier plot, but it's not actually necessary
plt.style.use('ggplot')

"""
This is to ignore RuntimeWarning: invalid value encountered in true_divide
I know that when my populations are zero there's some division by zero and
the resulting error terminates my function, which I want to avoid in this case.
"""
np.seterr(divide='ignore', invalid='ignore')

# These are the parameter values we'll be using
a = 0.005
b = 0.5
c = 0.5
d = 0.1
h = 0.1
K = 2000

# Define the system of ODEs
# P[0] is prey, P[1] is predator
def fish(P, t=0):
    return ([b*P[0]*(1-P[0]/K) - (a*P[0]*P[1])/(1+a*h*P[0]),
            c*(a*P[0]*P[1])/(1+a*h*P[0]) - d*P[1] ])

# Define equilibrium point
EQ = ([d/(a*(c-d*h)),b*(1+a*h*(d/(a*(c-d*h))))*(1-(d/(a*(c-d*h)))/K)/a])

"""
I need to define the possible values my initial points will take as they
relate to the equilibrium point. In this case I chose to plot 10 trajectories
ranging from 0.1 to 5
"""
values = np.linspace(0.1, 5, 10)
# I want each trajectory to have a different color
vcolors = plt.cm.autumn_r(np.linspace(0.1, 1, len(values)))

# Open figure
f = plt.figure()
"""
I need to define a range of time over which to integrate the system of ODEs
The values don't really matter in this case because our system doesn't have t
on the right hand side of dx/dt and dy/dt, but it is a necessary input for
integrate.odeint.
"""
t = np.linspace(0, 150, 1000)

# Plot trajectories by looping through the possible values
for v, col in zip(values, vcolors):
    # Starting point of each trajectory
    P0 = [E*v for E in EQ]
    # Integrate system of ODEs to get x and y values
    P = integrate.odeint(fish, P0, t)
    # Plot each trajectory
    plt.plot( P[:,0], P[:,1],
            # Different line width for different trajectories (optional)
            lw=0.5*v,
            # Different color for each trajectory
            color=col,
            # Assign starting point to trajectory label
            label='P0=(%.f, %.f)' % ( P0[0], P0[1]) )
"""
To plot the direction fields we first need to define a grid in order to
compute the direction at each point
"""
# Get limits of trajectory plot
ymax = plt.ylim(ymin=0)[1]
xmax = plt.xlim(xmin=0)[1]
# Define number of points
nb_points = 20
# Define x and y ranges
x = np.linspace(0, xmax, nb_points)
y = np.linspace(0, ymax, nb_points)
# Create meshgrid
X1 , Y1 = np.meshgrid(x,y)
# Calculate growth rate at each grid point
DX1, DY1 = fish([X1, Y1])
# Direction at each grid point is the hypotenuse of the prey direction and the
# predator direction.
M = (np.hypot(DX1, DY1))
# This is to avoid any divisions when normalizing
M[ M == 0] = 1.
# Normalize the length of each arrow (optional)
DX1 /= M
DY1 /= M

plt.title('Trajectories and direction fields')
"""
This is using the quiver function to plot the field of arrows using DX1 and
DY1 for direction and M for speed
"""
Q = plt.quiver(X1, Y1, DX1, DY1, M, pivot='mid', cmap=plt.cm.plasma)
plt.xlabel('Prey abundance')
plt.ylabel('Predator abundance')
plt.legend(bbox_to_anchor=(1.05, 1.0))
plt.grid()
plt.xlim(0, xmax)
plt.ylim(0, ymax)
plt.show()

This should produce the following plot. All P0s are the initial conditions we defined.
trajectories

We can also see that this parameter combination produces limit cycles in our system. If we change the parameter values to:

a = 0.005
b = 0.5
c = 0.5
d = 0.1
h = 0.1
K = 200

i.e. reduce the available resources to the prey, our trajectories look like this:

trajectories1

The equilibrium becomes stable, attracting the trajectories to it.

The same can be seen if we increase the predator death rate:

a = 0.005
b = 0.5
c = 0.5
d = 1.5
h = 0.1
K = 2000

trajectories2

The implication of this observation is that an initially stable system, can become unstable given more resources for the prey or less efficient predators. This has been referred to as the Paradox of Enrichment and other predator-prey models have tried to address it (more on this in future posts).

P.S: I would also like to link to this scipy tutorial, that I found very helpful and that contains more plotting tips.

Setting up Python and Eclipse

According to its website, Python is:

…an easy to learn, powerful programming language. It has efficient high-level data structures and a simple but effective approach to object-oriented programming. Python’s elegant syntax and dynamic typing, together with its interpreted nature, make it an ideal language for scripting and rapid application development in many areas on most platforms.

The Python interpreter and the extensive standard library are freely available in source or binary form for all major platforms from the Python Web site,http://www.python.org/, and may be freely distributed. The same site also contains distributions of and pointers to many free third party Python modules, programs and tools, and additional documentation.

This post covers how to set up Python and the Eclipse development environment.  We also provide a collection of posts on how to use Python for data analysis, starting here.

————————————————————————————————————————————————————–

PYTHON:

————————————————————————————————————————————————————–

The first step is to download Python and its various packages that will likely be useful to you at some point.

Python itself is available at: http://www.python.org/download/

I would recommend downloading and installing version 2.7.2, the latest production release under the 2.X series.  Also, stick with the 32-bit version as most all packages will be available for this version.  Avoid Python 3.X for now.  It is not as widely supported among the various Python packages that you might find useful and as such, should be avoided for now.  Keep in mind that there are some syntax differences as well between versions 2.X and 3.X that would need to be addressed whenever it does come time to update.

Just use the default settings during installation.

NOTE: If you have Cygwin installed on your system, it too has likely installed a version of Python.  Whenever you run Python from the command line, you should be careful to ensure that you are using the version that you expect (i.e., the default Cygwin installed Python versus the one that you installed).  Just be aware of this.  In general, it is easy to identify the version being picked up from the path name.  Also, it is generally best to use the version that you have installed.  It will usually be located in C:\Python27 whereas the Cygwin version will be located in C:\Cygwin\bin.

Now, install the various packages that may be useful. You should always be careful to install a version of the package that matches your version of Python (i.e., 2.7 if you are following my instructions).  Sometimes, if a package is not available for the version you are using (i) you may still be able to use it, or (ii) you may need to make minor tweaks to the package source to get things running. Also, always download the package installers, not the source.  Here are the common ones that you should definately install:

  • NumPy and SciPy available at http://numpy.scipy.org/.  These packages are useful for performing scientific computing within Python.  Download the “win32 superpacks” for each of these packages for the version of Python that you have installed.
  • PIL – the Python Imaging Library available at http://www.pythonware.com/products/pil/.  This package is useful to manipulating image files.
  • matplotlib – a 2D plotting library with Matlab-like syntax available at http://matplotlib.sourceforge.net/.  This package is very good for creating good publication quality figures.  If you starting using it, you will probably notice that the appearance of the figures, even on-screen, is much improved over what Matlab can produce.

The following are some optional packages based on your particular needs:

  • Py2exe – a package for bundling Python scripts into MS Windows executable programs available at http://www.py2exe.org/.  This is what I use to bundle all of the libraries and source code required by AeroVis into a self contained package that can be installed on any Windows system without the need to build or install Python, VTK, Qt, etc.
  • wxPython – GUI package for Python available at http://wxpython.org/.  Note, this is for developing graphical user interfaces (GUIs) for your Python scripts, it is not a GUI for Python.
  • PyQt – another GUI package for Python available at http://www.riverbankcomputing.co.uk/software/pyqt/intro.  PyQt is a set bindings for Nokia’s Qt application framework – a very rich and full featured graphical interface development framework.  AeroVis uses PyQt for its graphical interface.

————————————————————————————————————————————————————–

ECLIPSE:

————————————————————————————————————————————————————–

Now that you have Python and all of your needed packages installed, you can now move on to Eclipse. Eclipse is available from http://www.eclipse.org/downloads/packages/release/indigo/r.  The latest release (and probably the version you should be using) is Indigo.  Since we primarily use Visual Studio for C/C++ development, I would recommend downloading the IDE for Java as this will serve to provide you with a Java environment should you choose to explore this down the road.  I think you should be able to install either the 32-bit or 64-bit versions without issue.  Just make sure you are running a 64-bit OS if you choose to install that version.  When you go to download, Penn State actually has a mirror so choose this.  BTW, don’t choose the BitTorrent option – not a good idea on PSU networks.

Once you have downloaded the zip file containing Eclipse, you just unzip it wherever you want it to be installed.  This includes portable drives etc.  The beauty of Eclipse is that unlike many Windows programs, it is completely self contained and as such, can be run from any location.  Once unzipped, create a shortcut to the Eclipse executable and start it up.

————————————————————————————————————————————————————–

PYDEV:

————————————————————————————————————————————————————–

Now that Eclipse is installed, we can add a Python development environment inside Eclipse that will provide a very nice Python IDE with debugging capabilities, etc.

The install for packages inside Eclipse proceeds a little differently than what you may be used to.

The best option for installing PyDev is probably to install Aptana Studio which includes a variety of development tools.  Go to this site for instructions http://www.aptana.com/downloads/start or read on.

1) In the Eclipse Help menu, select Install New Software
2) Paste this URL into the Work With box: http://download.aptana.com/studio3/plugin/install
3) Check the box for Aptana Studio and click Next
4) Accept the license, etc., and restart Eclipse

Another option is to only install PyDev from within Eclipse, carefully follow the instructions available at: http://pydev.org/manual_101_install.html.  There’s no need for me to rehash all of these instructions here as they are quite good at the PyDev site.

Once PyDev is installed, you should be ready to go.

————————————————————————————————————————————————————–

Let me know if you run into any problems by leaving a comment.

————————————————————————————————————————————————————–

Up Next Time…

Developing and debugging Python scripts and projects in Eclipse