Predator-Prey Models and Global Warming Followup

I thought I'd follow up my last post with some information on some of the tools I used. I'm not going to talk about how I made the equations, because they came out terrible and I'm trying to come up with a new way. I guess if anyone is curious: I made them in LaTeX and used gimp to convert them to PNGs.

Ok, so let's see how I made the graphics. Here's the first graphic, in case anyone forgot.

The first step in creating the graphic was, obviously, to solve the ODEs. I used SciPy for that. Specifically, I used the odeint command. It's pretty simple to use in fact. Here's what it looks like:

	def derivative_fxn(x, time, ...):
	dx = ...
	return dx
	initial_values = array(...)
	times = arange(t0, tfinal, interval)
	odeint(derivative_fxn, initial_values, times, args=(...), ...)

view raw gistfile1.py hosted with ❤ by GitHub

This function is part of the NumPy python packages, so you can easily work with systems of equations, like the one seen in Predator-Prey Models and Global Warming. The function derivative_fxn needs to return the derivative given x and the time. Note, that x may be an array. The args=(...) allows you to pass a tuple of arguments to give to the derivative_fxn.

Anyway, so that is odeint. Now, for the full script used to solve the predator-prey models:

	#!/usr/bin/python
	# -*- coding: utf-8 -*-
	from scipy import *
	from scipy.integrate import odeint
	#Set the constants.
	D0 = 1.07173*10**12
	S0= 5.35866*10**10
	kb = 8.6173324*10**(-5)
	K = 1./21
	tlo = 20
	thi = 30
	#Formating string for printing the numbers
	formatStr = "%8.3f"
	#Death rate as a function of temperature
	def D(T):
	return D0 * exp(-0.66/(kb * (T + 273.15)))
	#Sexy time as a function of temperature
	def S(T):
	return S0 * exp(-0.66/(kb * (T + 273.15)))
	#Calculate the derivative
	def derivative(species, time, temp):
	crow_der = S(temp) * species[0] - K * species[0] * species[1]
	cat_der = K * species[0] * species[1] - D(temp) * species[1]
	return (crow_der, cat_der)
	#Set the time intervals and initial values
	t = arange(0,16,1*10**(-2))
	species_i = array([110,7])
	#Go through the temperatures
	temps = range(tlo, thi + 1)
	print temps
	result = empty( (len(t),2,len(temps)) )
	#For each temperature, solve the ODE
	for i,temp in zip(range(len(temps)), temps):
	result[:,:,i] = odeint(derivative, species_i, t, args=(temp,), printmessg=True)
	#Print them out
	lines = []
	for i in range(len(t)):
	lines.append(formatStr % t[i])
	for j in range(len(temps)):
	lines.append(" ")
	lines.append(formatStr % result[i,0,j])
	for j in range(len(temps)):
	lines.append(" ")
	lines.append(formatStr % result[i,1,j])
	lines.append("\n")
	with open("out.txt", "w") as f:
	f.write("".join(lines))

view raw gistfile1.py hosted with ❤ by GitHub

The advantage to using SciPy is the ability to spice up the ODE with some programming, which allowed me to calculate the solution at each of the temperatures pretty easily. Now, how to graph it? I've been pretty hot on R lately for my graphs. Some people like the combination of SciPy and matplotlib, which is powerful. But, for publication quality graphics, R is what I use.

I won't go into too much detail on the script, but you can probably figure out most of it from below. The great thing about python and R is that you can open up the languages in an interactive terminal and go line by line to examine how a script works. The disadvantage to being an interactive language, especially for R, is the performance hit. Anyway, without further ado here is the R plotting script:

	data <- read.table("out.txt")
	time <- data[,1]
	temps <- 20:30
	crows <- data[,2:(length(temps) + 1)]
	cats <- data[,(length(temps) + 2):ncol(data)]
	colors <- colorRampPalette(c("blue", "red")) (length(temps))
	cairo_pdf("crows.pdf", width=5, height=7)
	par(family="LMRoman10", fg="dark gray", mfrow=c(2,1), mar=c(2,4,0.5,2))
	plot(time, crows[,1], ylim=c(min(crows), max(crows)), xlab="", ylab="Crows", lwd=1, col=colors[1], type="l", xaxt="n")
	for(i in 2:length(temps)){
	lines(time, crows[,i], col=colors[i], ylim=c(min(crows), max(crows)))
	}
	par(mar=c(5,4,0,2))
	plot(time, cats[,1], ylim=c(min(cats), max(cats)), xlab="time", ylab="Cats", lwd=1, col=colors[1], type="l")
	for(i in 2:length(temps)){
	lines(time, cats[,i], col=colors[i], ylim=c(min(cats), max(cats)))
	}
	graphics.off()

view raw gistfile1.r hosted with ❤ by GitHub

That's it! Now I just need to figure out how to format equations...