Using SAGE to investigate the discrete logistic equation
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)