[Tutor] Lotka-Volterra Model Simulation Questions
Oscar Benjamin
oscar.j.benjamin at gmail.com
Sat Sep 29 00:46:08 CEST 2012
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>
More information about the Tutor
mailing list