from vpython import gcurve, color
# line 17, foxRate formula was incorrect
# line 27, expression for foxNew referenced rab instead of fox
# line 32, initial point needs to be a tuple with time at 0
# line 33, initial point needs to be a tuple with time at 0
# line 38, population is second tuple element, not first
# line 39, population is second tuple element, not first
# line 41, needs to be indented to be part of the loop
# line 44, needs to be indented to be part of the loop
# line 45, needs to be indented to be part of the loop
# line 51, fox population should be plotted in blue
# line 52, argument should be foxPoints, not rabPoints
# parameters
alpha = 0.1 # intrinsic rate of prey population increase
beta = 0.01 # predation rate coefficient
gamma = 0.05 # predator mortality rate
delta = 0.001 # reproduction rate of predators per 1 prey eaten
rab0 = 50 # initial rabbit population
fox0 = 15 # initial predator population
timelimit = 400 # number of timesteps to graph
# Lotka-Volterra equations
def rabRate(rab, fox):
return rab * (alpha - beta * fox)
def foxRate(rab, fox):
return -fox * (gamma - delta * rab)
# takes current populations and returns tuple of new populations
def pop(rab, fox):
# second-order Runge-Kutta method
# ensure estimated population is never negative
rabMid = max(rab + 0.5 * rabRate(rab, fox), 0)
foxMid = max(fox + 0.5 * foxRate(rab, fox), 0)
rabNew = max(rab + rabRate(rabMid, foxMid), 0)
foxNew = max(fox + foxRate(rabMid, foxMid), 0)
return rabNew, foxNew
# initialize time-population lists with the initial populations at time 0
rabPoints = [(0, rab0)]
foxPoints = [(0, fox0)]
# calculate new populations and add the data onto the lists
for t in range(1, timelimit):
# get current populations
rabCur = rabPoints[-1][1] # last list element, 2nd element
foxCur = foxPoints[-1][1] # last list element, 2nd element
# calculate new populations
rabNext, foxNext = pop(rabCur, foxCur)
# lists are of tuples with time (x-axis) followed by population (y-axis)
rabPoints.append((t, rabNext))
foxPoints.append((t, foxNext))
# plot populations versus time
# rabbits in orange, foxes in blue
rabPlot = gcurve(color=color.orange)
rabPlot.plot(pos=rabPoints)
foxPlot = gcurve(color=color.blue)
foxPlot.plot(pos=foxPoints)