[Tutor] Lotka-Volterra Model Simulation Questions

[Oscar Benjamin]

On 28 September 2012 21:32, Jim Apto <jimmyapt0 at gmail.com> wrote:

> Hello folks,
>
> I'm relatively new to python, and was asked to program a lotka-volterra
> model (predator and prey relation) simulator.  The program basically will
> basically have a menu that takes user input, collect data, and then create
> a graph.  Currently i've been working on the simulator section; I can't
> seem to get the lists right.  I've assigned the following variables and
> parameters to the model for the program:
>
> x represents prey population
> y represents predator population
> dy/dt and dx/dt represents growth rate of the two populations over time
> t represents time
>
> a is the growth rate of prey
> b is the rate at which predators kill prey
> g is the death rate of predators
> d is the rate at which the predators population increases by consuming prey
>
> The equation:
> dx/dt = x(a-by)
> dy/dt = -y(g-dx)
>
> The code I have for this section is:
> def deltaX(a,b,x,y):
>     dx = x*(a-b*y)
>
> def deltaY(g,d,x,y):
>     dy = -y*(g-d*x)
>

The normal way to program an ODE solver is to define a function that takes
a vector input X and returns a vector dX/dt. The idea is that rather than
keeping a separate state for x and y you keep a combined state [x, y]. This
makes sense since the state of your system at any time is defined by both
values. Keeping that in mind I would write a function like

def derivative(Z):
    x, y = Z
   dxdt = ...
   dydt = ...
   dZdt = [dxdt, dydt]
   return dZdt

This function uses lists of numbers to represent the state vector. You can
then plug this into a function that uses a particular ODE solving algorithm:

def euler(f, Z1, dt):
    dZdt = f(Z1)
    Z2 = []
    for i in range(len(Z1)):
        Z2.append(Z1[i] + dt * dZdt[i])
    return Z2

Or a shorter version:

def euler(f, Z1, dt):
    return [z + dt * dz for z, dz in zip(Z, f(Z))]

You then get the new state by plugging the old state into the euler
function repeatedly, e.g.:

dt = .001
Z0 = [0.5, 0.5]   # Initial condition
Z1 = euler(derivative, Z0, dt)
Z2 = euler(derivative, Z1, dt)
...

Once you can get all the simulation values you will be in a position to
think about how to rearrange the lists so that you can plot them.

By the way, this sort of thing is usually done with numpy (that makes a few
of these things a bit easier).

Oscar
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://mail.python.org/pipermail/tutor/attachments/20120928/dce617c6/attachment.html>