# -*- coding: utf-8 -*- # 3.0 # # Population with several species - Systems of ODE. # #from sympy import * #from sympy.utilities.lambdify import lambdify #init_printing() import numpy as np from scipy.integrate import odeint import matplotlib.pyplot as plt import math # # A linear model for two species. # #(b1, b2) = symbols('b1 b2') #(c1, c2) = symbols('c1 c2') #(d1, d2) = symbols('d1 d2') #t = Symbol('t') #p = Function('p')(t) #q = Function('q')(t) #f1 = b1+c1*p+d1*q #f2 = b2+c2*p+d2*q #twolin = [Eq(p.diff(t),f1),Eq(q.diff(t),f2)] #pprint(twolin) # sympy dsolve doesn't solve systems of ODEs, but it will be implemented in the next release (0.7.6) def vectorfield(w, t, c): """ Defines the differential equations for the population dynamics. Arguments: w : vector of the species variables: w = [p,q] t : time c : vector of the parameters: c = [b1,b2,c1,c2,d1,d2] """ p, q = w b1, b2, c1, c2, d1, d2 = c # Create f = (p',q'): f = [b1+c1*p+d1*q, b2+c2*p+d2*q] return f # # Example: arm race - steady scenario. # # Parameter values b1 = 1.5; b2 = 0.0 c1 = -1.0; c2 = 0.5 d1 = 0.5; d2 = -1.0 # Initial conditions p0 = 0.5 q0 = 1.5 # ODE solver parameters abserr = 1.0e-8 relerr = 1.0e-6 stoptime = 10.0 numpoints = 250 t = [stoptime * float(i) / (numpoints - 1) for i in range(numpoints)] # Pack up the parameters and initial conditions: c = [b1,b2,c1,c2,d1,d2] w0 = [p0,q0] # Call the ODE solver. linsol = odeint(vectorfield, w0, t, args=(c,), atol=abserr, rtol=relerr) # plot the solution p, q = linsol.T plt.figure() plt.plot(t,p,t,q) plt.show() # # Example: arms race - unlimited growth. # # Parameter values b1 = 0.5; b2 = -5./4. c1 = -3./4.; c2 = 1.0 d1 = 1.0; d2 = -3./4. # Initial conditions p0 = 1.0 q0 = 2.5 # Pack up the parameters and initial conditions: c = [b1,b2,c1,c2,d1,d2] w0 = [p0,q0] # Call the ODE solver. linsol = odeint(vectorfield, w0, t, args=(c,), atol=abserr, rtol=relerr) # plot the solution p, q = linsol.T plt.figure() plt.plot(t,p,t,q) plt.show() # # Example: the peaceful neighbour # # Parameter c2 is negative, so country 2 will destroy weapons, if its neighbour's armament grows. # # Parameter values b1 = 0.0; b2 = 5./2. c1 = -3./4.; c2 = -1.0 d1 = 1.0; d2 = -3./4. # Initial conditions p0 = 0.5 q0 = 3.0 # Pack up the parameters and initial conditions: c = [b1,b2,c1,c2,d1,d2] w0 = [p0,q0] # Call the ODE solver. linsol1 = odeint(vectorfield, w0, t, args=(c,), atol=abserr, rtol=relerr) # plot the solution p, q = linsol1.T plt.figure() plt.plot(t,p,t,q) plt.show() # # A non-linear model for two species. # def vectorfield_nonlin(w, t, c): """ Defines the differential equations for the population dynamics. Arguments: w : vector of the species variables: w = [p,q] t : time c : vector of the parameters: c = [b1,b2,c1,c2,d1,d2] """ p, q = w b1, b2, c1, c2, d1, d2 = c # Create f = (p',q'): f = [(b1+c1*p+d1*q)*p, (b2+c2*p+d2*q)*q] return f # # Example: Competition - Steady State. # # Competition between two species (sharing common ressources, for example). # # Parameter values b1 = 2.5 + math.sqrt(3.)/24.; b2 = 7./8. + math.sqrt(27.)/2. c1 = -5./8.; c2 = -math.sqrt(27.)/8. d1 = -math.sqrt(3.)/24.; d2 = -7./8. # Initial conditions p0 = 0.25 q0 = 3.0 # ODE solver parameters abserr = 1.0e-8 relerr = 1.0e-6 stoptime = 10.0 numpoints = 250 t = [stoptime * float(i) / (numpoints - 1) for i in range(numpoints)] # Pack up the parameters and initial conditions: c = [b1,b2,c1,c2,d1,d2] w0 = [p0,q0] # Call the ODE solver. nonlinsol = odeint(vectorfield_nonlin, w0, t, args=(c,), atol=abserr, rtol=relerr) # plot the solution p, q = nonlinsol.T plt.figure() plt.plot(t,p,t,q) plt.show() # # Example: Competition - Extinction. # # In this example, one species dies out. # # Parameter values b1 = 71./8.; b2 = 73./8. c1 = -23./12.; c2 = -25./12. d1 = -25./12.; d2 = -23./12. # Initial conditions p0 = 0.25 q0 = 0.5 # Pack up the parameters and initial conditions: c = [b1,b2,c1,c2,d1,d2] w0 = [p0,q0] # Call the ODE solver. nonlinsol = odeint(vectorfield_nonlin, w0, t, args=(c,), atol=abserr, rtol=relerr) # plot the solution p, q = nonlinsol.T plt.figure() plt.plot(t,p,t,q) plt.show() # # A non-linear model for 2 species (predator-prey): Volterra-Lotka. # # Parameter values b1 = -0.5; b2 = 1./5. c1 = 0.0; c2 = -1./50. d1 = 1./200.; d2 = 0.0 # Initial conditions p0 = 6.0 q0 = 50.0 # Pack up the parameters and initial conditions: c = [b1,b2,c1,c2,d1,d2] w0 = [p0,q0] stoptime = 100.0 numpoints = 500 t = [stoptime * float(i) / (numpoints - 1) for i in range(numpoints)] # Call the ODE solver. nonlinsol = odeint(vectorfield_nonlin, w0, t, args=(c,), atol=abserr, rtol=relerr) # plot the solution p, q = nonlinsol.T plt.figure() plt.plot(t,p,t,q) plt.show()