Qualitative analysis and Bifurcation diagram Tutorial
Posted: | Source
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:
