# -*- coding: utf-8 -*- # 3.0 # # Model of Maltus # from sympy import * from sympy.utilities.lambdify import lambdify from sympy.core.relational import * init_printing() # # Define the model of Maltus as a differential equation in sympy. # alpha = Symbol('alpha') t = Symbol('t') p = Function('p')(t) maltusODE = Eq(p.diff(t),alpha*p) pprint(maltusODE) # # Let sympy compute the general solution (note: the solution is given as an equation!) # maltusSol = dsolve(maltusODE,p).rhs pprint(maltusSol) # C1 = Symbol('C1') maltusFun = maltusSol.subs([(C1,1),(alpha,-0.5)]) pprint(maltusFun) # import numpy as np from mpl_toolkits.mplot3d import Axes3D import matplotlib.pyplot as plt from matplotlib import cm lam_maltus = lambdify((t), maltusFun) Nt = 100 T = np.linspace(0.,10.,Nt) P = [lam_maltus(Te) for Te in T] plt.figure() plot(T, P) plt.show() # # .. and an example of a growing population: # maltusFun = maltusSol.subs([(C1,-0.5),(alpha,0.3)]) pprint(maltusFun) # lam_maltus = lambdify((t), maltusFun) Nt = 100 T = np.linspace(0.,10.,Nt) P = [lam_maltus(Te) for Te in T] plt.figure() plot(T, P) plt.show() # # Model of Verhulst # beta = Symbol('beta') verhulstODE = Eq(p.diff(t),alpha-beta*p) pprint(verhulstODE) # verhulstSol = dsolve(verhulstODE,p).rhs pprint(verhulstSol) # # First example: the saturation value is alpha/beta = 2, and the initial population # is larger than the saturation value: # verhulstFun = simplify(verhulstSol.subs([(C1,ln(0.5)),(alpha,1.0),(beta,0.5)])) pprint(verhulstFun) # lam_verhulst = lambdify((t), verhulstFun) Nt = 100 T = np.linspace(0.,10.,Nt) P = [lam_verhulst(Te) for Te in T] plt.figure() plot(T, P) plt.show() # # Second example: initial population is smaller than saturation value: # verhulstFun = simplify(verhulstSol.subs([(C1,ln(-0.5)),(alpha,1.0),(beta,0.5)])) pprint(verhulstFun) # lam_verhulst = lambdify((t), verhulstFun) Nt = 100 T = np.linspace(0.,10.,Nt) P = [lam_verhulst(Te) for Te in T] plt.figure() plot(T, P) plt.show() # # Logistic Growth # beta = Symbol('beta') logisticODE = Eq(p.diff(t),alpha*(1-p/beta)*p) pprint(logisticODE) # logisticSol = dsolve(logisticODE,p,hint='best').rhs pprint(logisticSol) # logisticFun = simplify(logisticSol.subs([(C1,9.5),(alpha,1.0),(beta,2.0)])) pprint(logisticFun) # lam_logistic = lambdify((t), logisticFun) Nt = 100 T = np.linspace(0.,10.,Nt) P = [lam_logistic(Te) for Te in T] plt.figure() plot(T, P) plt.show() # logisticFun = simplify(logisticSol.subs([(C1,-1./6.),(alpha,1.0),(beta,2.0)])) pprint(logisticFun) # lam_logistic = lambdify((t), logisticFun) Nt = 100 T = np.linspace(0.,10.,Nt) P = [lam_logistic(Te) for Te in T] plt.figure() plot(T, P) plt.show() # # Logistic Growth with Threshold # delta = Symbol('delta') thresholdODE = Eq(p.diff(t),alpha*(1-p/beta)*(1-p/delta)*p) pprint(thresholdODE) # thresholdSol = dsolve(thresholdODE,p,hint='all') print(thresholdSol) # # sympy, at the moment of writing this document, does not provide the solution of the logistic growth with thresholds as a symbolic expersion, but rather as an equation, which we cannot use directly for plotting. Therefore, we will solve the ODE numerically to demonstrate results. # #thresholdODE_1 = thresholdODE.subs([(alpha,-1),(beta,2),(delta,4)]) #pprint(thresholdODE_1) #thresholdSol = dsolve(thresholdODE_1,p,hint='best') #pprint(thresholdSol) # def threshold_dwdt(w, t, p): """ Compute time derivative for the logistic growth with thresholds model. Arguments: w : species variable t : time p : vector of the parameters: p = [alpha,beta,delta] """ a, b, d = p # Create f = w': f = a*(1-w/b)*(1-w/d)*w return f # import numpy as np from scipy.integrate import odeint # symbolic solution with specified coefficients thresholdODE_1 = thresholdODE.subs([(alpha,-1),(beta,2),(delta,4)]) pprint(thresholdODE_1) thresholdSol = dsolve(thresholdODE_1,p,hint='best') pprint(thresholdSol) # Parameter values a = -1. b = 2. d = 4. # Initial conditions p0 = 2.2 # 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 = [a,b,d] # Call the ODE solver. linsol = odeint(threshold_dwdt, p0, T, args=(c,), atol=abserr, rtol=relerr) # plot the solution P = linsol plt.figure() plt.plot(T,P) plt.show() # # Initial conditions p0 = 1.9 # Call the ODE solver. linsol = odeint(threshold_dwdt, p0, T, args=(c,), atol=abserr, rtol=relerr) # plot the solution P = linsol plt.figure() plt.plot(T,P) plt.show() # # Initial conditions p0 = 5.0 # Call the ODE solver. linsol = odeint(threshold_dwdt, p0, T, args=(c,), atol=abserr, rtol=relerr) # plot the solution P = linsol plt.figure() plt.plot(T,P) plt.show()