"""
We assume we have two species, herbivores with population x, 
and predators with propulation y. We assume that x grows 
exponentially in the absence of predators, and that y decays 
exponentially in the absence of prey. Consider, say, the system

x' = x - xy/2, y' = -y + xy/2.          (*)

We call this the predator-prey model or the Lotka-Volterra model.

The following Sage commands implement this.

REFERENCES:
[1] J. Rosenberg's Matlab page, 
    http://www.math.umd.edu/~jmr/246/predprey.html
[2] Wikipedia http://en.wikipedia.org/wiki/Lotka-Volterra_model

"""

#This plots the well-known "phase portrait" 
#  {(x,y) | x = x(t), y = y(t) satisfy (*)}.
#  (x0,y0)= (0.5,0.5) is in red
#  (x0,y0)= (0.75,0.75) is in gree
#  (x0,y0)= (1.0, 1.0) is in blue
# See predator-prey1.png
t, x, y = PolynomialRing(QQ,3,"txy").gens()
f = x - x*y/2
g = -y + x*y/2
pts = eulers_method_2x2(f,g, 0.0, 0.5, 0.5, 0.01, 8.0, method=None)
pts_xy = [[x[1],x[2]] for x in pts]
P1 = line(pts_xy,rgbcolor=(1,0,0))
pts = eulers_method_2x2(f,g, 0.0, 0.75, 0.75, 0.01, 8.0, method=None)
pts_xy = [[x[1],x[2]] for x in pts]
P2 = line(pts_xy,rgbcolor=(0,1,0))
pts = eulers_method_2x2(f,g, 0.0, 1.0, 1.0, 0.01, 8.0, method=None)
pts_xy = [[x[1],x[2]] for x in pts]
P3 = line(pts_xy,rgbcolor=(0,0,1))
show(P1+P2+P3)

#This plots the periodic function x = x(t) which satisfies (*).
#  (x0,y0)= (0.5,0.5) is in red
#  (x0,y0)= (0.75,0.75) is in gree
#  (x0,y0)= (1.0, 1.0) is in blue
# See predator-prey2.png
pts = eulers_method_2x2(f,g, 0.0, 0.5, 0.5, 0.001, 20.0, method=None)
pts_tx = [[x[0],x[1]] for x in pts]
P1 = line(pts_tx,rgbcolor=(1,0,0))
pts = eulers_method_2x2(f,g, 0.0, 0.75, 0.75, 0.001, 20.0, method=None)
pts_tx = [[x[0],x[1]] for x in pts]
P2 = line(pts_tx,rgbcolor=(0,1,0))
pts = eulers_method_2x2(f,g, 0.0, 1.0, 1.0, 0.001, 20.0, method=None)
pts_tx = [[x[0],x[1]] for x in pts]
P3 = line(pts_tx,rgbcolor=(0,0,1))
show(P1+P2+P3)