#!/usr/bin/env python
# coding: utf-8

# > This is one of the 100 recipes of the [IPython Cookbook](http://ipython-books.github.io/), the definitive guide to high-performance scientific computing and data science in Python.
# 

# # 15.7. Analyzing a nonlinear differential system: Lotka-Volterra (predator-prey) equations

# Here, we conduct a brief analytical study of a famous nonlinear differential system: the Lotka-Volterra equations, also known as predator-prey equations. This simple model describes the evolution of two interacting populations (e.g. sharks and sardines), where the predators eat the preys. This example illustrates how we can use SymPy to obtain exact expressions and results for fixed points and their stability.

# In[ ]:


from sympy import *
init_printing()


# In[ ]:


var('x y')
var('a b c d', positive=True)


# The variables x and y represent the populations of the preys and predators, respectively. The parameters a, b, c and d are positive parameters (described more precisely in "How it works..."). The equations are:
# 
# $$\begin{align}
# \frac{dx}{dt} &= f(x) = x(a-by)\\
# \frac{dy}{dt} &= g(x) = -y(c-dx)
# \end{align}$$

# In[ ]:


f = x * (a - b*y)
g = -y * (c - d*x)


# Let's find the fixed points of the system (solving f(x,y) = g(x,y) = 0).

# In[ ]:


solve([f, g], (x, y))


# In[ ]:


(x0, y0), (x1, y1) = _


# Let's write the 2D vector with the two equations.

# In[ ]:


M = Matrix((f, g)); M


# Now we can compute the Jacobian of the system, as a function of (x, y).

# In[ ]:


J = M.jacobian((x, y)); J


# Let's study the stability of the two fixed points by looking at the eigenvalues of the Jacobian at these points.

# In[ ]:


M0 = J.subs(x, x0).subs(y, y0); M0


# In[ ]:


M0.eigenvals()


# The parameters a and c are strictly positive, so the eigenvalues are real and of opposite signs, and this fixed point is a saddle point. Since this point is unstable, the extinction of both populations is unlikely in this model.

# In[ ]:


M1 = J.subs(x, x1).subs(y, y1); M1


# In[ ]:


M1.eigenvals()


# The eigenvalues are purely imaginary so this fixed point is not hyperbolic, and we cannot draw conclusions about the qualitative behavior of the system around this fixed point from this linear analysis. However, one can show with other methods that oscillations occur around this point.

# > You'll find all the explanations, figures, references, and much more in the book (to be released later this summer).
# 
# > [IPython Cookbook](http://ipython-books.github.io/), by [Cyrille Rossant](http://cyrille.rossant.net), Packt Publishing, 2014 (500 pages).