[SciPy-user] nullclines for nonlinear ODEs

[John Hunter]

In a workshop on scientific computing in python  I ran with Andrew
Straw last weekend at the Claremont Colleges, I presented a unit on
ODEs, numerically integrating the Lotka-Volterra equations, drawing a
direction field, plotting a trajectory in the  phase plane, and
plotting the nullclines.  For the nullclines, I used matplotlib's
contouring routine with levels=[0], but I was wondering if there was a
better/more direct route.

The example script is included below (and attached)

import numpy as n
import pylab as p
import scipy.integrate as integrate

def dr(r, f):
    return alpha*r - beta*r*f

def df(r, f):
    return gamma*r*f - delta*f

def derivs(state, t):
    """
    Map the state variable [rabbits, foxes] to the derivitives
    [deltar, deltaf] at time t
    """
    #print t, state
    r, f = state  # rabbits and foxes
    deltar = dr(r, f)  # change in rabbits
    deltaf = df(r, f) # change in foxes
    return deltar, deltaf

alpha, delta = 1, .25
beta, gamma = .2, .05

# the initial population of rabbits and foxes
r0 = 20
f0 = 10

t = n.arange(0.0, 100, 0.1)

y0 = [r0, f0]  # the initial [rabbits, foxes] state vector
y = integrate.odeint(derivs, y0, t)
r = y[:,0]  # extract the rabbits vector
f = y[:,1]  # extract the foxes vector

p.figure()
p.plot(t, r, label='rabbits')
p.plot(t, f, label='foxes')
p.xlabel('time (years)')
p.ylabel('population')
p.title('population trajectories')
p.grid()
p.legend()
p.savefig('lotka_volterra.png', dpi=150)
p.savefig('lotka_volterra.eps')


p.figure()
p.plot(r, f)
p.xlabel('rabbits')
p.ylabel('foxes')
p.title('phase plane')


# make a direction field plot with quiver
rmax = 1.1 * r.max()
fmax = 1.1 * f.max()
R, F = n.meshgrid(n.arange(-1, rmax), n.arange(-1, fmax))
dR = dr(R, F)
dF = df(R, F)
p.quiver(R, F, dR, dF)


R, F = n.meshgrid(n.arange(-1, rmax, .1), n.arange(-1, fmax, .1))
dR = dr(R, F)
dF = df(R, F)

p.contour(R, F, dR, levels=[0], linewidths=3, colors='black')
p.contour(R, F, dF, levels=[0], linewidths=3, colors='black')
p.ylabel('foxes')
p.title('trajectory, direction field and null clines')

p.savefig('lotka_volterra_pplane.png', dpi=150)
p.savefig('lotka_volterra_pplane.eps')


p.show()
-------------- next part --------------
A non-text attachment was scrubbed...
Name: lotka_volterra.py
Type: text/x-python
Size: 1706 bytes
Desc: not available
URL: <http://mail.scipy.org/pipermail/scipy-user/attachments/20071030/cb74a054/attachment.py>