Log in

No account? Create an account

Previous Entry | Next Entry

As you may know I am a biology teacher. Next term I will teach Ecology. One part of that is the relationship between predators and prey (wolves and rabbits, for example). There are quite some mathematical model for these relationships, in school we mostly teach a simplified version of the Lotka-Volterra system.

This kind of stuff really needs a simple application to test the influence of certain factors. For example,

  • what happens if after two years we kill all but 5 predators?

  • what happens if after two years we kill all but 5 prey?

  • what happens if after 1 year, for whatever reasons, the rabbits become more fertile?

There are tools schools can buy, but a) they suck and b) they are very expensive.

Introducing: "Predator and Prey" by Carsten Niehaus (tm)

Free Image Hosting at www.ImageShack.us


As you can see the basic stuff all works. When happens when I start with just on predator but 20 prey but let the predators be more efficient hunters? Look here:

Free Image Hosting at www.ImageShack.us


The application is written in PyQt4 using matplotlib. I would really like to express one thing: Python + Qt4 rocks. Seriously, the whole application has only 306 lines including comments!

Are you interested in the code or even in helping me writing a cool PP-Simulation? Do you think this should go into KDE EDU? I am hosting everything in GitHub: PP in GitHub.


( 11 comments — Leave a comment )
Jul. 2nd, 2009 02:38 pm (UTC)
I just changed the plot style
The lines were not thin enough, here (and with 20000 iterations (the other pics show 5000)) your are:

Free Image Hosting at www.ImageShack.us

QuickPost Quickpost this image to Myspace, Digg, Facebook, and others!
Jul. 2nd, 2009 04:44 pm (UTC)
After looking at the screenshots I think there is something wrong with your simulation. Both populations of an Lotka-Volterra system should be cyclic while yours increase in time.

After having a quick look at the sources, I can even explain why:
The integrator you use (the method to integrate the Differential equation over time) is not appropriate for the problem.
The method you are using is known as Euler Method (http://en.wikipedia.org/wiki/Euler_method), which is quite straightforward but sadly also quite bad because it produces huge errors (your populations seem to increase exponentially). A more appropriate Integrator would be for example a Runge-Kutta-Integrator of second order.

I have an implementation of several integrators in python but they use the numpy package which you probably dont want as dependency. But if you want I can rewrite one of them for your purpose.
I will send you an email with an patch if you are interested.
Jul. 2nd, 2009 05:14 pm (UTC)
Re: Integrator
I am very much interested, yes!

cniehaus kde org
Jul. 2nd, 2009 05:17 pm (UTC)
Re: Integrator
I just noticed that Alpha even knows those equations! What do you think about:

Jul. 2nd, 2009 06:55 pm (UTC)
Re: Integrator
The formulas given by wolfram alpha are actually already in your program (implemented in the functions PredatorPreyCalculator.dx and PredatorPreyCalculator.dy). These are the differential equations that define the Lotka-Volterra System.
The difficult part is to calculate the actual trajectories (the long-term behaviour) from these. As the computer simulates the system using finite time-steps, there will always be errors. These can be minimized by using very small time-steps, but they will never vanish (In your program the size of the timestep is determined by PredatorPreyCalculator.dt and the errors will become smaller the smaller you set this value). But using small timesteps also means that simulating the system takes longer, therefor it is necessary to search for other means to reduce the errors. Such a method is using an integrator that performs better for the same size of timesteps.

The easiest integrator is the euler integrator which is a first order Taylor expansion:
"x_new = x_old + stepsize * dx(x_old,T)"
where dx(x,T) is the differential equation defining the system (aka the derivative of x)

The problem with this approach can be easily seen when you try to calculate the exponential function with it (using the differential equation "dx(x) = x"):
For example with a stepsize of 1 and an x_old=1=exp(0):
"exp(1) = x_new = 1 + 1 * 1 = 2"
Where you would expect x_new = exp(1) = e (eulers number) =~ 2.7
The problem can be lessened by setting the stepsize to 0.5. But the you will need two steps to calculate the value of exp(1) as you have to calculate exp(0.5) first. And still: the value you gain for exp(1) will be to small.

For this reason one usually uses more sophisticated integrators today. One of the still very simple ones is a Runge-Kutta of second order. I have sent you the modified "PredatorPreyCalculator.py" using this integration method per mail, I hope it works.

Sorry for the way too long comment.

Jul. 3rd, 2009 08:06 am (UTC)
Re: Integrator
Wow, thanks for that explanation. I never thought about all this (well, perhaps the last time at university a couple of years ago ;-).

Thanks for your code, it works very well!
Jul. 3rd, 2009 08:58 am (UTC)
Re: Integrator
I was looking at http://docs.scipy.org/doc/scipy/reference/integrate.html and indeed do not see runga kutta 4rth order, only lsoda which has adaptive time stepping and error control and often a better choice than rk4.

I'm writing an extension on the ode integrators in scipy to have more of the functionality as found in matlab/maple/...:

I'd be glad to add some educational codes like euler, or rk4. If you have a pure python implementation, it is bsd licenced like python/scipy, then send it to me 'Stephan'.

Showing how code must be used in very important, I added the doublependulum with output an mpeg movie in the docs:
An example of the predator prey would be nice too. :-) As this is educational, can you licence it bsd? I am all for gpl, but when it comes to these kind of things, bsd or public domain is the norm.
Jul. 3rd, 2009 12:38 pm (UTC)
Re: Integrator
look here: http://github.com/cniehaus/biotools/blob/c84a67065f524f60851a574cf08f7c482af0442a/raeuber-beute/predatorpreycalculator.py

That is pure Python. I will relicence it to BSD, just need to check for the exact licencetext.

Jul. 3rd, 2009 06:46 pm (UTC)
Re: Integrator
Done, you will find the relicenced file in github.
Jul. 3rd, 2009 02:57 am (UTC)
Cool - I never thought of using QT4 with python before, but next time I need to prototype a GUI I think I'll take some notes from your code - thanks!

Oh, and PP sims rock, it was one of the first things that interested me about artificial life.
Jul. 3rd, 2009 08:06 am (UTC)
I will try to document the code very well. I not only teach biology and chemistry but also computer sciences. So a well documented code can as well be used as a demo for certain things.

For example, look at this commit:


Stephans code works way better, but a math teacher might want to demonstrate why the Euler methode is not the perfect choice here. This way he would just need to uncomment and comment a couple of lines and students could direclty SEE why RK2 is better here.
( 11 comments — Leave a comment )