COMP 130 Elements of Algorithms and Computation
Spring 2012

Even More Refining of the Predator-Prey Algorithm

The following is the Comp130 continuation of the "Improving the Accuracy of the Algorithm" notes from Comp200.

Lists of tuples

In the code from class, we created three lists:  times, prey (hare) population and predator (lynx) population.   While this may seem like an easy and "natural" solution, there are some issues to think about:

  1. Each population value, prey or predator, is intimately associated with a specific time value.
  2. Population values have no meaning without their corresponding time sequence.
  3. In the end, there are no guarantees that the number of elements in all three lists is the same.
  4. In the end, there are no guarantees that any given population value correctly associates with the time value with the same list index.
  5. The fact that the plotting functions want the data in a particular format should have no bearing on the output format of the populations function.

In effect, the initmately coupled nature of the triplet of (time, prey  population, predator population) that, while created together, was lost because the values were split into separate lists.

Solution:  change the loop to create a single list of tuples, keeping the associated time and population values together.    When needed, use the zip() function to transpose the list of tuples (or lists) into a list of lists for plotting.

Class Activity 1:  Implement this improvement.

Encapsulate the loop body

Looking at our code, we see that all the main loop is doing is to walk through the time steps, generating the next hare and lynx populations.   Notice that the looping process itself does not depend on

Also the functional process to calculate the next population values does not depend on

In such, we should encapsulate the body of the loop such that it simply and directly expresses what the loop does: repetitively generates the next population value.  

# global variables to define the various rate constants
hareBirthRate = .4 
harePredationRate = .003
lynxBirthRate = .004
lynxDeathRate = .2 	

def populations(hares0, lynxes0, nYears, stepsPerYear):
	"""
	Generate a list of tuples, (t, hares, lynxes), starting at the initial 
	populations, hares0 and lynxes0 for nYears at stepsPerYears times per year.   
	The initial time, t=0 and thetime nYear are both included in the output.
	"""
	popl = (0.0, hares0, lynxes0)   # initial population at time = 0

	results = [popl]    # initialize the results list
	dt = 1.0/stepsPerYear  # calculate the time step
	for i in range(nYears*stepsPerYear):    # loop over all steps
		popl = nextPopl(popl, dt)  # calculate the next population given the current population
		results.append(popl)  # append to the results
		print popl
	return results
	

Note that the rate constants are not passed into the nextPopl function as parameters.   Is this good or bad?

Pros:

Cons:

Is there a way of getting the pros but not cons?   The problem is that we don't yet have the techniques to properly encapsulate the rate constants with the nextPopl function.   Stay tuned!

Click  here for advanced Plotting How-To's  (Direct link to download plotfuncs.py

Class Activity 2:

  1. Write the the nextPopl() function.
    def nextPopl(popl, dt):
    	""" Calculate the next population based on the current (given) population and a time step, 
    	using a generalized coupled differential equation.
    	dHares/dt == Psi(hares, lynxes)
    	dLynxes/dt == Phi(hares, lynxes)
    	where
    	popl == (time, hares, linxes) == the current hare and lynx populations at specified time
    	dt == time step in units of years 
    	Returns a tuple of the next population values. 
    	"""
    	# calculate the first derivative
    	dHare_dt = psi(popl)
    	dLynx_dt = phi(popl)
    	
    	return (popl[0]+dt,  popl[1] + dHare_dt*dt , popl[2] + dLynx_dt*dt ) 
    
    
    def psi(popl):
    	""" Calculate the net hare population growth rate
    	popl is a tuple of the current populations: (time, hares, lynxes)
    	dHares/dt == Psi(hares, lynxes)
    	""" 
    	return popl[1] * (hareBirthRate - harePredationRate*popl[2])
    	
    def phi(popl):
    	""" Calculate the net lynx population growth rate
    	popl is a tuple of the current populations: (time, hares, lynxes)
    	dLynxes/dt == Phi(hares, lynxes)
    	""" 
    	return popl[2]*(lynxBirthRate*popl[1]- lynxDeathRate)		
    		
  2. Implement the whole algorithm using  your nextPopl() function, the above looping code and a plotting function.
    import plotfuncs as pf  # this should be at top of file.  Need to have plotfuncs.py in the same folder.
    
    initHares = 100
    initLynxes = 50
    nYears = 50
    stepsPerYear = 12
    		
    popData1 = populations(initHares, initLynxes, nYears, stepsPerYear)
    pf.makeTimePlot(popData1, (('Hares', 'o'), ('Lynxes', 's')), 'Years', 'Population', 'Predator-Prey Populations')	
    

Class Activity 3:

  1. By simply rearranging your results from Class Activity 2, and using an x-y plotting function, create a graph of the predator population vs. the prey-population.
    results = zip(*popData1)   # transpose the data into separate time, hare and lynx population lists
    xydata = []   # list of x-y plot data lists
    xydata.append(zip(results[1], results[2]))   # this can be done multiple times to add more x-y graph data lists.
    pf.makeXYPlot(xyData, (('prey vs. predator dataset #1', '')), 'Hare Population', 'Lynx Population', 'Predator-Prey Populations')	
    		
  2. Create graphs to show how the predator-prey population cycles change as various parameters are changed.

 

2nd Order Approximations

The problem with our predator-prey algorithm is that is uses a linear ("first-order") approximation to find the next population value given the current population value:

P(t+dt) ≈ P(t) + slope * dt  = P(t) + (dP/dt) * dt

But what if P(t) is not a straight line, that is, dP/dt is not a constant?   If this is true, then the above approximation must fundamentally give the wrong results.  Luckily, as we saw, if the time step, dt, is small enough, the population errors can be made arbitrarily small, though at the cost of increased computation time and resources.

If the slope of the population curve changes in time, then for any finite time interval,  dt, the actual population value will be different than the first-order result we calculated:

2nd order approx.

 

But back-tracking a bit, the above first-order formula is actually part of a more general approximation scheme from calculus called a "Taylor expansion" or "Taylor series" or "Maclaurin series":   

Our population approximation is really just the first two terms of the Taylor expansion, the "zero'th order" term, which is simply the value of P(t), the original point, and the "first order" term, which is the slope (derivative) of P(t) with respect to time multiplied by the time difference to the next point.    But this ignores all the rest of the terms of expansion, which account for how the slope changes in time and how the slope of the slope changes in times, and so on.  In the limit of an infinite number of terms being used, the Taylor expansion becomes an exact result.

We can approximate the curvature of the P(t)  line by simply adding in the next term of the Taylor expansion, the second order term which is the slope of the slope, i.e. the second derivative of the population with respect to time:

P(t+dt) ≈  P(t) + (dP/dt) * dt  + ½ (d2P/dt2) dt2

Those that have taken elementary physics should recognize the resulting equation from the equations of motion of an object subject to a constant acceleration (change of velocity per time):

x(t+dt) =  x(t) + v * dt  + ½ a dt2

where v = dx/dt and a = d(dx/dt)/dt = d2x/dt2.

It is important to realize that a Taylor expansion enables you to appoximate the value of a function at some later time, t+dt, given only values that are calculated at the present time, t.

Continue on with the 2nd-order Lotka-Volterra model.