py33. Population with Harvesting in Python

A harvest term, H, is used in the equation of the population, which is:

dN/dt = rN(1-N/K) - H

Steady-state solutions can be found from the roots of equation above, after setting dN/dt equal to zero. We find the roots using numpy.roots function.

# ex33.py
from __future__ import print_function, division
import numpy as np
import matplotlib.pyplot as plt

r = .10 # growth rate / yr
K = 100000 # carrying capacity
t = 100 # number of years
H = 1600 # Harvesting / yr
num = np.zeros(t+1)
num0 = np.zeros(t+1)
num[0] = 40000
num0[0] = 40000
for i in range(t):
    num[i+1] = num[i]+r*num[i]*(1-num[i]/K)-H
    num0[i+1] = num0[i]+r*num0[i]*(1-num0[i]/K)
sol = np.roots([r/K,-r,H])
print('\nSolutions are %.0fk and %.0fk.' %
      (sol[0]/1000,sol[1]/1000))
mStr = ("Growth rate: %.2f,"
        "Carrying Capacity: %dk"
        "\nHarvesting %d")
tStr = r,K/1000,H
main =  mStr % tStr
plt.plot(range(t+1),num/1000, 'b',
         range(t+1),num0/1000,'m')
plt.ylim(38,102)
plt.xlabel('Year')
plt.ylabel('Number (Thousands)')
plt.title(main)
plt.text(70, 84, '80K (H = %d)' % H)
plt.text(70, 97, '100K (H = 0)')

#Solutions are 80k and 20k.

Output: