## Using SAGE to investigate the discrete logistic equation

December 10th, 2009 | Categories: Fractals, math software, Open Source, sage interactions | Tags:

There have been a couple of blog posts recently that have focused on creating interactive demonstrations for the discrete logistic equation – a highly simplified model of population growth that is often used as an example of how chaotic solutions can arise in simple systems.  The first blog post over at Division by Zero used the free GeoGebra package to create the demonstration and a follow up post over at MathRecreation used a proprietary package called Fathom (which I have to confess I have never heard of – drop me a line if you have used it and like it).  Finally, there is also a Mathematica demonstration of the discrete logistic equation over at the Wolfram Demonstrations project.

I figured that one more demonstration wouldn’t hurt so I coded it up in SAGE – a free open source mathematical package that has a level of power on par with Mathematica or MATLAB.  Here’s the code (click here if you’d prefer to download it as a file).

def newpop(m,prevpop):
return m*prevpop*(1-prevpop)

def populationhistory(startpop,m,length):
history = [startpop]
for i in range(length):
history.append( newpop(m,history[i]) )
return history

@interact
def _( m=slider(0.05,5,0.05,default=1.75,label='Malthus Factor') ):
myplot=list_plot( populationhistory(0.1,m,20) ,plotjoined=True,marker='o',ymin=0,ymax=1)
myplot.show()

Here’s a screenshot for a Malthus Factor of 1.75

and here’s one for a Malthus Factor of 3.1

There are now so many different ways to easily make interactive mathematical demonstrations that there really is no excuse not to use them.

Update (11th December 2009)

As Harald Schilly points out, you can make a nice fractal out of this by plotting the limit points.  My blogging software ruined his code when he placed it in the comments section so I reproduce it here (but he’s also uploaded it to sagenb)

var('x malthus')
step(x,malthus) = malthus * x * (1-x)
stepfast = fast_callable(step, vars=[x, malthus], domain=RDF)

def logistic(m, step):
# filter cycles
v = .5
for i in range(100):
v = stepfast(v,m)
points = []
for i in range(100):
v = stepfast(v,m)
points.append((m+step*random(),v))
return points
points=[]
step = 0.005
for m in sxrange(2.5,4,step):
points += logistic(m, step)
point(points,pointsize=1).show(dpi=150)

1. If you just plot the limit points you get a nice fractal ;)

cell 1:

var(‘x malthus’)
step(x,malthus) = malthus * x * (1-x)
stepfast = fast_callable(step, vars=[x, malthus], domain=RDF)

def logistic(m, step):
# filter cycles
v = .5
for i in range(100):
v = stepfast(v,m)
points = []
for i in range(100):
v = stepfast(v,m)
# points with horizontal jitter
points.append((m+step*random(),v))
return points

cell 2:

points=[]
step = 0.005
for m in sxrange(2.5,4,step):
points += logistic(m, step)
point(points,pointsize=1).show(dpi=150)

2. Ok, that’s ugly here… published on sagenb.org
http://sagenb.org/pub/1240/

3. You might be interested in these interacts if you haven’t seen them already:

http://wiki.sagemath.org/interact/dynsys

Cheers,
Marshall

4. Hi Marshall

So your logistic map interact is ‘cythonized’ – pretend I am an idiot for a while (not too much of a stretch I know)……

Cythonized means that those functions run as compiled-C right? Or put another way…very VERY fast?

Cheers,
Mike

5. Hi Mike,
Thanks for blogging about this. I’ve done something very similar as a lab in Calc I (not Cythonized, though, so probably slower than Marshall’s), and for pedagogical use I found it helpful to make sliders for the starting population and the step size as well. This particularly shows (as with other Euler-type approximations) that the ‘chaotic’ coefficients and behavior depend a lot on the step size, but not so much on the starting point (with certain exceptions). Of course, given your description in ‘about’, probably you won’t need it too much for teaching :) but anyway it’s fun to think about how to use these capabilities.
Best,
– kcrisman

6. I’ve made this (the limit-y, fractal-y) demo a million times in Mathematica a million different ways, and every time I see it I love it more :D

7. @Connor – a million different ways huh? Sounds interesting! Care to post a couple? ;)

@Kcrisman – I don’t teach maths but I do teach math software – often to lecturers or postgrads. I like examples such as this because they show what can be achieved with very little actual effort