R33. Population with Harvesting in R

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 to zero. We use the function polyroot to find the roots.

# ex33.R
r <- .10 # growth rate / yr
K <- 100000 # carrying capacity
t <- 100 # number of years
H <- 1600 # Harvesting / yr
num <- num0 <- numeric(t+1)
num[1] <- num0[1] <- 40000
for (i in 1: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 <- polyroot(c(H,-r,r/K))
cat('\nSolutions are ',
    Re(sol[1])/1000, 'k and ',
    Re(sol[2])/1000, 'k.\n',
    sep = '')
main <- paste("Growth rate: ",r,
              ", Carrying Capacity: ",
              K/1000,"k\nHarvesting: ", H, sep = '')
plot(num/1000, type = 'l', col = 'blue',
     xlab = 'Year', ylab = 'Number (Thousands)',
     main = main, ylim =c(num[1]/1000,K/1000))
lines(num0/1000, col = 'magenta' )
text(x = 80, y = 84, paste('80K (H =',H,')'), col = 'blue')
text(x = 80, y = 97, '100K (H = 0)', col = 'magenta')

#Solutions are 20k and 80k.

Output: