Import packages
#Importation des packages
import numpy as np
from scipy.integrate import odeint
import matplotlib.pyplot as plt
Table 1
#Données prédateurs
time_data = [0,5,10,15,20,25,30,35,40,50]
predator_data = [265,245,200,65,8,0,0,0,0,0]
Table 2
#Données proies
deer_density_data = [0,0.005,0.01,0.015,0.02,0.025,0.05]
predation_rate_data = [0,3,13,28,51,56,56]
Table 3
#Données des ressources
food_per_deer_data = [0,500,1000,1500,2000,10000]
growth_rate_data = [-0.5,-0.15,0,0.15,0.2,0.2]
Table 4
#Données de végétation
veg_density_data = [0,0.25,0.5,0.75,1]
reg_time_data = [35,15,5,1.5,1]
Initial values (BC) & Time parameters
# Initial conditions (account for boundaries counditions for differential eqt°)
D0 = 4000 # initial deer pop
F0 = 4.7*(10**8) # initial food
Z0 = [D0, F0] # list of initial values
# Simulation parameters (time)
ti = 0 # initial time
tf = 50 # final time
nt = 1001 # number of time points
t = np.linspace(ti,tf,nt) # vector of time points
Model function
# function that returns dZdt = [dD/dt ; dF/dt]
def model(Z,t):
# parameters
area = 800000
max_food_capacity = 4.8*(10**8)
daily_requirement = 2000
feeding_period = 1
# subsetting du vecteur input Z
D = Z[0]
F = Z[1]
# Deer dynamics
deer_density = D / area
predation_rate = np.interp(deer_density, deer_density_data, predation_rate_data)
predator = np.interp(t, time_data, predator_data)
predation_loss = predation_rate*predator
food_per_deer = F/D
growth_rate = np.interp(food_per_deer,food_per_deer_data, growth_rate_data)
net_increase = D*growth_rate
# Food dynamics
food_demand = daily_requirement * D
if food_demand < F/feeding_period :
browsing_loss = food_demand
else :
browsing_loss = F/feeding_period
veg_density = F / max_food_capacity
regeneration_time = np.interp(veg_density, veg_density_data, reg_time_data)
food_increase = (max_food_capacity - F)/regeneration_time
# computation des dérivées
dDdt = net_increase - predation_loss
dFdt = food_increase - browsing_loss
dZdt = [dDdt,dFdt]
return dZdt
# Time evolution
Z = odeint(model,Z0,t)
D = Z[:,0]
F = Z[:,1]
Visualisation (plot)
# plot timeseries
plt.plot(t,D,'--b',label='Deer')
plt.xlabel('time [years]')
plt.legend(loc='best')
ax.spines['right'].set_visible(False)
ax.spines['top'].set_visible(False)
plt.show()
Execution error
NameError: name 'ax' is not defined
plt.plot(t,F,'r',label='Food')
ax = plt.gca()
plt.xlabel('time [years]')
plt.legend(loc='best')
ax.spines['right'].set_visible(False)
ax.spines['top'].set_visible(False)
plt.show()
(Les deux en un ... )
fig, ax = plt.subplots(1, 1)
ax.plot(t,D,'--b',label='Deer')
ax2 = ax.twinx()
ax2.plot(t,F,'r',label='Food')
plt.draw()
State trajectory
plt.plot(F,D)
plt.xlabel('food [kcal*year/day]')
plt.ylabel("Number of deers")
plt.legend(loc='best')
ax.spines['right'].set_visible(False)
ax.spines['top'].set_visible(False)
plt.show()
No handles with labels found to put in legend.
Without predation loss (Q2b)
# function that returns dZdt = [dD/dt ; dF/dt]
def model(Z,t):
# parameters
area = 800000
max_food_capacity = 4.8*(10**8)
daily_requirement = 2000
feeding_period = 1
# subsetting du vecteur input Z
D = Z[0]
F = Z[1]
# Deer dynamics
deer_density = D / area
food_per_deer = F/D
growth_rate = np.interp(food_per_deer,food_per_deer_data, growth_rate_data)
net_increase = D*growth_rate
# Food dynamics
food_demand = daily_requirement * D
if food_demand < F/feeding_period :
browsing_loss = food_demand
else :
browsing_loss = F/feeding_period
veg_density = F / max_food_capacity
regeneration_time = np.interp(veg_density, veg_density_data, reg_time_data)
food_increase = (max_food_capacity - F)/regeneration_time
# computation des dérivées
dDdt = net_increase
dFdt = food_increase - browsing_loss
dZdt = [dDdt,dFdt]
return dZdt
# Time evolution
Z = odeint(model,Z0,t)
D = Z[:,0]
F = Z[:,1]
fig, ax = plt.subplots(1, 1)
ax.plot(t,D,'--b',label='Deer')
ax2 = ax.twinx()
ax2.plot(t,F,'r',label='Food')
plt.draw()