Qualitative analysis and Bifurcation diagram Tutorial
Posted:
Qualitative analysis and Bifurcation diagram Tutorial
This tutorial assumes you have read the tutorial on numerical integration.
Exploring the parameter space: bifurcation diagrams
The Rosenzweig-MacArthur consumer-resource model
\[ \begin{aligned} \frac{dR}{dt} &= rR \left( 1 - \frac{R}{K} \right) - \frac{a R C}{1+ahR} \\ \frac{dC}{dt} &= \frac{e a R C}{1+ahR} - d C \end{aligned} \]
Rosenzweig–MacArthur model solutions
%matplotlib inline
from numpy import *
from scipy.integrate import odeint
from matplotlib.pyplot import *
ion()
def RM(y, t, r, K, a, h, e, d):
return array([ y[0] * ( r*(1-y[0]/K) - a*y[1]/(1+a*h*y[0]) ),
y[1] * (e*a*y[0]/(1+a*h*y[0]) - d) ])
t = arange(0, 1000, .1)
y0 = [1, 1.]
pars = (1., 10., 1., 0.1, 0.1, 0.1)
y = odeint(RM, y0, t, pars)
plot(t, y)
xlabel('time')
ylabel('population')
legend(['resource', 'consumer'])
The phase space flow and the fixed point
# plot the solution in the phase space
plot(y[:,0], y[:,1])
# defines a grid of points
R, C = meshgrid(arange(0.95, 1.25, .05), arange(0.95, 1.04, 0.01))
# calculates the value of the derivative at the point in the grid
dy = RM(array([R, C]), 0, *pars)
# plots arrows on the points of the grid, with the difection
# and length determined by the derivative dy
# This is a picture of the flow of the solution in the phase space
quiver(R, C, dy[0,:], dy[1,:], scale_units='xy', angles='xy')
xlabel('Resource')
ylabel('Consumer')
Messing a little with the parameters...
# now K = 15
t = arange(0, 1000, 1.)
pars = (1., 15., 1., 0.1, 0.1, 0.1)
y_osc = odeint(RM, y0, t, pars)
plot(t, y_osc)
xlabel('time')
ylabel('population')
legend(['resource', 'consumer'])
plot(y_osc[:,0], y_osc[:,1])
R, C = meshgrid(arange(0, 6., .4), arange(0, 2.1, 0.2))
dy = RM(array([R, C]), 0, *pars)
quiver(R, C, dy[0,:], dy[1,:], scale_units='xy', angles='xy')
xlabel('R')
ylabel('C')
plot(10., y[-500:,0].min(), 'og')
plot(10., y[-500:,0].max(), 'og')
plot(10., y[-500:,1].min(), 'ob')
plot(10., y[-500:,1].max(), 'ob')
plot(15., y_osc[-500:,0].min(), 'og')
plot(15., y_osc[-500:,0].max(), 'og')
plot(15., y_osc[-500:,1].min(), 'ob')
plot(15., y_osc[-500:,1].max(), 'ob')
xlim((0, 20))
yscale('log')
xlabel('K')
ylabel('min / max population')
The paradox of enrichment
What happens when we change the carrying capacity \(K\) from very small values up to very large values? For very small values, the resource is not going to sustain the consumer population, but for larger values ok \(K\), both species should be benefited... right?
# empty lists to append the values later
ymin = []
ymax = []
KK = arange(.5, 25, .5)
t = arange(0, 6000, 1.)
# loop over the values of K (KK)
for K in KK:
# redefine the parameters using the new K
pars = (1., K, 1., 0.1, 0.1, 0.1)
# integrate again the equation, with new parameters
y = odeint(RM, y0, t, pars)
# calculate the minimum and maximum of the populations, but
# only for the last 1000 steps (the long-term solution),
# appending the result to the list
# question: is 1000 enough? When it wouldn't be?
ymin.append(y[-1000:,:].min(axis=0))
ymax.append(y[-1000:,:].max(axis=0))
# convert the lists into arrays
ymin = array(ymin)
ymax = array(ymax)
# and now, we plot the bifurcation diagram
plot(KK, ymin[:,0], 'g')
plot(KK, ymax[:,0], 'g')
plot(KK, ymin[:,1], 'b')
plot(KK, ymax[:,1], 'b')
xlabel('$K$')
ylabel('min/max populations')
legend(['resource', 'consumer'])
# use a log scale in the y-axis
yscale('log')
Well, the first prediction was OK (notice that the plot above uses a log scale), but for high \(K\), the minima of the oscillation go to very low values, so that the populations have a high risk of extinction.
Consumer-resource dynamics in a seasonal environment
\[ \begin{aligned} \frac{dR}{dt} &= r(t) R \left( 1 - \frac{R}{K} \right) - \frac{a R C}{1+ahR} \\ \frac{dC}{dt} &= \frac{e a R C}{1+ahR} - d C \\ r(t) &= r_0 (1+\alpha \sin(2\pi t/T)) \end{aligned} \]
def RM_season(y, t, r, alpha, T, K, a, h, e, d):
# in this function, `t` appears explicitly
return array([ y[0] * ( r * (1+alpha*sin(2*pi*t/T)) *
(1-y[0]/K) - a*y[1]/(1+a*h*y[0]) ),
y[1] * (e*a*y[0]/(1+a*h*y[0]) - d) ])
t = arange(0, 2000, 1.)
y0 = [1., 1.]
pars = (1., 0.1, 80., 10., 1., 0.1, 0.1, 0.1)
y = odeint(RM_season, y0, t, pars)
plot(t, y)
xlabel('time')
ylabel('population')
legend(['resource', 'consumer'])
A resonance diagram
ymin = []
ymax = []
t = arange(0, 6000, 1.) # times
TT = arange(1, 80, 2) # periods
for T in TT:
pars = (1., 0.1, T, 10., 1., 0.1, 0.1, 0.1)
y = odeint(RM_season, y0, t, pars)
ymin.append(y[-1000:,:].min(axis=0))
ymax.append(y[-1000:,:].max(axis=0))
ymin = array(ymin)
ymax = array(ymax)
plot(TT, ymin[:,0], 'g')
plot(TT, ymax[:,0], 'g')
plot(TT, ymin[:,1], 'b')
plot(TT, ymax[:,1], 'b')
xlabel('$T$')
ylabel('min/max populations')
legend(['resource', 'consumer'])
yscale('log')
What if I really want to explore a 10-dimensional parameter space?
First: good luck. Second, you will probably have to sample the space, rather than go through the whole thing. The recommended method to do it is using so-called Latin Hypercube samples, that uses a random sampling while ensuring a roughly regularly-spaced distribution. There are implementations for both R and python:
Numerical Integration Tutorial
Posted:
Numerically solving differential equations with python
This is a brief description of what numerical integration is and a practical tutorial on how to do it in Python.
Software required
In order to run this notebook in your own computer, you need to install the following software:
- python
- numpy and scipy - python scientific libraries
- matplotlib - a library for plotting
- the ipython notebook
On Windows and Mac, we recommend installing the Anaconda distribution, which includes all of the above in a single package (among several other libraries), available at http://continuum.io/downloads.
On Linux, you can install everything using your distribution's prefered way, e.g.:
- Debian/Ubuntu:
sudo apt-get install python-numpy python-scipy python-matplotlib python-ipython-notebook - Fedora:
sudo yum install python-numpy python-scipy python-matplotlib python-ipython-notebook - Arch:
sudo pacman -S python-numpy python-scipy python-matplotlib ipython python-tornado python-jinja
Code snippets shown here can also be copied into a pure text file with .py extension and ran outside the notebook (e.g., in an python or ipython shell).
From the web
Alternatively, you can use a service that runs notebooks on the cloud, e.g. SageMathCloud or wakari. It is possible to visualize publicly-available notebooks using http://nbviewer.ipython.org, but no computation can be performed (it just shows saved pre-calculated results).
How numerical integration works
Let's say we have a differential equation that we don't know how (or don't want) to derive its (analytical) solution. We can still find out what the solutions are through numerical integration. So, how dows that work?
The idea is to approximate the solution at successive small time intervals, extrapolating the value of the derivative over each interval. For example, let's take the differential equation
\[ \frac{dx}{dt} = f(x) = x (1 - x) \]
with an initial value \(x_0 = 0.1\) at an initial time \(t=0\) (that is, \(x(0) = 0.1\)). At \(t=0\), the derivative \(\frac{dx}{dt}\) values \(f(0.1) = 0.1 \times (1-0.1) = 0.09\). We pick a small interval step, say, \(\Delta t = 0.5\), and assume that that value of the derivative is a good approximation over the whole interval from \(t=0\) up to \(t=0.5\). This means that in this time \(x\) is going to increase by \(\frac{dx}{dt} \times \Delta t = 0.09 \times 0.5 = 0.045\). So our approximate solution for \(x\) at \(t=0.5\) is \(x(0) + 0.045 = 0.145\). We can then use this value of \(x(0.5)\) to calculate the next point in time, \(t=1\). We calculate the derivative at each step, multiply by the time step and add to the previous value of the solution, as in the table below:
| \(t\) | \(x\) | \(\frac{dx}{dt}\) |
|---|---|---|
| 0 | 0.1 | 0.09 |
| 0.5 | 0.145 | 0.123975 |
| 1.0 | 0.206987 | 0.164144 |
| 1.5 | 0.289059 | 0.205504 |
| 2.0 | 0.391811 | 0.238295 |
Of course, this is terribly tedious to do by hand, so we can write a simple program to do it and plot the solution. Below we compare it to the known analytical solution of this differential equation (the logistic equation). Don't worry about the code just yet: there are better and simpler ways to do it!
%matplotlib inline
from numpy import *
from matplotlib.pyplot import *
# time intervals
tt = arange(0, 10, 0.5)
# initial condition
xx = [0.1]
def f(x):
return x * (1.-x)
# loop over time
for t in tt[1:]:
xx.append(xx[-1] + 0.5 * f(xx[-1]))
# plotting
plot(tt, xx, '.-')
ta = arange(0, 10, 0.01)
plot(ta, 0.1 * exp(ta)/(1+0.1*(exp(ta)-1.)))
xlabel('t')
ylabel('x')
legend(['approximation', 'analytical solution'], loc='best',)
Why use scientific libraries?
The method we just used above is called the Euler method, and is the simplest one available. The problem is that, although it works reasonably well for the differential equation above, in many cases it doesn't perform very well. There are many ways to improve it: in fact, there are many books entirely dedicated to this. Although many math or physics students do learn how to implement more sophisticated methods, the topic is really deep. Luckily, we can rely on the expertise of lots of people to come up with good algorithms that work well in most situations.
Then, how... ?
We are going to demonstrate how to use scientific libraries to integrate differential equations. Although the specific commands depend on the software, the general procedure is usually the same:
- define the derivative function (the right hand side of the differential equation)
- choose a time step or a sequence of times where you want the solution
- provide the parameters and the initial condition
- pass the function, time sequence, parameters and initial conditions to a computer routine that runs the integration.
A single equation
So, let's start with the same equation as above, the logistic equation, now with any parameters for growth rate and carrying capacity:
\[ \frac{dx}{dt} = f(x) = r x \left(1 - \frac{x}{K} \right) \]
with \(r=2\), \(K=10\) and \(x(0) = 0.1\). We show how to integrate it using python below, introducing key language syntax as necessary.
# everything after a '#' is a comment
## we begin importing libraries we are going to use
# import all (*) functions from numpy library, eg array, arange etc.
from numpy import *
# import all (*) interactive plotting functions, eg plot, xlabel etc.
from matplotlib.pyplot import *
# import the numerical integrator we will use, odeint()
from scipy.integrate import odeint
# time steps: an array of values starting from 0 going up to (but
# excluding) 10, in steps of 0.01
t = arange(0, 10., 0.01)
# parameters
r = 2.
K = 10.
# initial condition
x0 = 0.1
# let's define the right-hand side of the differential equation
# It must be a function of the dependent variable (x) and of the
# time (t), even if time does not appear explicitly
# this is how you define a function:
def f(x, t, r, K):
# in python, there are no curling braces '{}' to start or
# end a function, nor any special keyword: the block is defined
# by leading spaces (usually 4)
# arithmetic is done the same as in other languages: + - * /
return r*x*(1-x/K)
# call the function that performs the integration
# the order of the arguments is as below: the derivative function,
# the initial condition, the points where we want the solution, and
# a list of parameters
x = odeint(f, x0, t, (r, K))
# plot the solution
plot(t, x)
xlabel('t') # define label of x-axis
ylabel('x') # and of y-axis
# plot analytical solution
# notice that `t` is an array: when you do any arithmetical operation
# with an array, it is the same as doing it for each element
plot(t, K * x0 * exp(r*t)/(K+x0*(exp(r*t)-1.)))
legend(['approximation', 'analytical solution'], loc='best') # draw legend
We get a much better approximation now, the two curves superimpose each other!
Now, what if we wanted to integrate a system of differential equations? Let's take the Lotka-Volterra equations:
\[ \begin{aligned} \frac{dV}{dt} &= r V - c V P\\ \frac{dP}{dt} &= ec V P - dP \end{aligned}\]
In this case, the variable is no longer a number, but an array [V, P]. We do the same as before, but now x is going to be an array:
# we didn't need to do this again: if the cell above was run already,
# the libraries are imported, but we repeat it here for convenience
from numpy import *
from matplotlib.pyplot import *
from scipy.integrate import odeint
t = arange(0, 50., 0.01)
# parameters
r = 2.
c = 0.5
e = 0.1
d = 1.
# initial condition: this is an array now!
x0 = array([1., 3.])
# the function still receives only `x`, but it will be an array, not a number
def LV(x, t, r, c, e, d):
# in python, arrays are numbered from 0, so the first element
# is x[0], the second is x[1]. The square brackets `[ ]` define a
# list, that is converted to an array using the function `array()`.
# Notice that the first entry corresponds to dV/dt and the second to dP/dt
return array([ r*x[0] - c * x[0] * x[1],
e * c * x[0] * x[1] - d * x[1] ])
# call the function that performs the integration
# the order of the arguments is as below: the derivative function,
# the initial condition, the points where we want the solution, and
# a list of parameters
x = odeint(LV, x0, t, (r, c, e, d))
VV = arange(0.1, 1., 0.01)
for
# Now `x` is a 2-dimension array of size 5000 x 2 (5000 time steps by 2
# variables). We can check it like this:
print('shape of x:', x.shape)
# plot the solution
plot(t, x)
xlabel('t') # define label of x-axis
ylabel('populations') # and of y-axis
legend(['V', 'P'], loc='upper right')
An interesting thing to do here is take a look at the phase space, that is, plot only the dependent variables, without respect to time:
# `x[0,0]` is the first value (1st line, 1st column), `x[0,1]` is the value of
# the 1st line, 2nd column, which corresponds to the value of P at the initial
# time. We plot just this point first to know where we started:
plot(x[0,0], x[0,1], 'o')
print('Initial condition:', x[0])
# `x[0]` or (equivalently) x[0,:] is the first line, and `x[:,0]` is the first
# column. Notice the colon `:` stands for all the values of that axis. We are
# going to plot the second column (P) against the first (V):
plot(x[:,0], x[:,1])
plot(x2[:,0], x2[:,1])
xlabel('V')
ylabel('P')
Congratulations: you are now ready to integrate any system of differential equations! (We hope generalizing the above to more than 2 equations won't be very challenging).
For more info:
- Python tutorial (chapters 3 to 5 are specially useful).
- An introduction to Numpy
- Another one, covering a little bit more ground.
- The matplotlib gallery: all kinds of plots, with sample code to use.
From Consumer-Resource to single-species dynamics
Posted:
From consumer-resource models to single-species dynamics
We are interested in building simple consumer-resource models and exploring how single-species models emerge as we simplify the dynamics of the resource. We assume a simple Lotka-Volterra type equation, with a carrying capacity for the resource:
