Limit Cycles
Limit Cycles: the Thin Line Between Life and Unlife [Rachael Mansbach, LANL]
Hands-On
- Limit Cycles
- Belousov–Zhabotinsky reaction
lifeandunlife.pptxBelousov-Zhabotinsky Reaction.mp4trajectory (Converted).movBelousov–Zhabotinsky reaction
#1) Solve the Lotka-Volterra Equations#Question: What happens when you change the initial conditions #Initial conditions are the # rabbits, # foxes#What if there are no rabbits and no foxes?#No rabbits?#No foxes?#Why?#What about changing the timespan?
def LVeqn(t,P):# # Lotka-Volterra model equations# ----------# Parameters# ----------# t: "time"# P: "population" or concentration (vector)# # -------# Returns# -------# dXdt: numpy array, essentially how the concentrations # change with time# #initial parameters alpha = 2.0 beta = 1.0 delta = 1.0 gamma = 1.0 rabbits = P[0] foxes = P[1] drabbit = alpha * rabbits - beta * rabbits * foxes dfoxes = delta * rabbits * foxes - gamma * foxes dPdt = np.array([drabbit,dfoxes]) return dPdt
y0 = np.array([5,5]) #this is the initial conditionstmax = 60tspan = (0,tmax)sol = solve_ivp(LVeqn,tspan,y0,method='BDF')ys = sol.yts = sol.t
#plot figure 1print("plotting fig 1")plt.figure()plt.plot(ts,ys[0,:])plt.plot(ts,ys[1,:],'--')plt.xlabel('time')plt.ylabel('population size')plt.legend(['rabbits','foxes'])plt.show()
#plot figure 2#comment out below this line (*) and above (**) to see Fig 1fig = plt.figure()ax = fig.add_subplot(111)print("plotting fig 2, comment out to see fig 1.")for i in range(0,tmax,5): ax.plot(ys[0,i:(i+10)],ys[1,i:(i+10)]) ax.plot(ys[0,i],ys[1,i],'r.') plt.annotate(str(i),(ys[0,i],ys[1,i])) #plt.pause(0.1)plt.xlabel('rabbits')plt.ylabel('foxes')plt.show()
#(**)
##2) Solve the Belousov-Zhabotinsky equations# These are a more complicated set of equations but they have similar# behavior# What happens when you change the initial conditions?# What happens if you change the parameters in the equation?# (try f = 0.5, eps = 0.2)# (This corresponds to changing concentrations of other chemicals in the mixture)
def BZeqn(t,P):# # BZ model equations# # ----------# Parameters# ----------# t: "time"# P: "population" or concentration (vector)# # -------# Returns# -------# dXdt: numpy array, essentially how the concentrations # change with time# x = P[0] z = P[1] #set of parameters in equation (external conditions) f = 1.0 eps = 0.2 q = 8e-4 #model equations dzdt = x-z dxdt = (1/eps)*(x*(1-x) + f*((q-x)/(q+x))*z) dXdt = np.array([dxdt,dzdt]) return dXdt
tmax = 200
y0 = np.array([0.2,0.5]) #initial conditionstspan = (0,tmax)sol = solve_ivp(BZeqn,tspan,y0,method='BDF')ys = sol.yts = sol.t
#plot figure 1print("plotting Fig 1")plt.figure()plt.plot(ts,ys[0,:])plt.plot(ts,ys[1,:],'--')plt.show() #plot figure 2, comment out below here (*) and above (**) to see Fig 1fig = plt.figure()ax = fig.add_subplot(1, 1, 1)print("plotting Fig 2; comment out to see Fig 1")for i in range(0,tmax,5): plt.plot(ys[0,i:(i+10)],ys[1,i:(i+10)]) plt.plot(ys[0,i],ys[1,i],'r.') plt.annotate(str(i),(ys[0,i],ys[1,i]))
#(**)
#Bonus: Competitive Lotka-Volterra Equations# What happens if you change the initial conditions just a tiny bit?
def compLVeqn(t,P):# # Very special set of "competitive" Lotka-Volterra# equations. Instead of predator-prey, this set of equations describes# two species that compete for the same resource, like two different# species of rabbits who both eat the same grass, although# in this particular case, we're working with a set of four species# dPdt = np.zeros(len(P)) r = np.array([1,0.72,1.53,1.27]) a = np.array([[1,1.09,1.52,0], [0,1,0.44,1.36], [2.33,0,1,0.47], [1.21,0.51,0.35,1]]) K = np.array([1,1,1,1]) for i in range(len(P)): s = 0 for j in range(len(P)): s += a[i,j] * P[j] dPdt[i] = r[i] * P[i] * (1 - s / K[i]) return dPdt
tmax = 200
y0 = np.array([10,8,2,.1])#initial conditionstspan = (0,tmax)sol = solve_ivp(compLVeqn,tspan,y0,method='BDF')ys = sol.yts = sol.ttmax = len(ts)
#figure 1print("plotting fig 1")plt.figure()plt.plot(ts,ys[0,:])plt.plot(ts,ys[1,:],'--')plt.plot(ts,ys[2,:],'-.')plt.plot(ts,ys[3,:],':')plt.show() #figure 2#comment out below this line (*) and above (**) to see fig 1fig = plt.figure()ax = fig.add_subplot(1, 1, 1)print("plotting fig 2; comment out to see fig 1")for i in range(0,tmax,5): plt.plot(ys[0,i:(i+10)],ys[1,i:(i+10)]) plt.plot(ys[0,i],ys[1,i],'r.') plt.annotate(str(i),(ys[0,i],ys[1,i]))
plt.show()#(**)