Instructions

Complete the following exercises either by hand or using RStudio. All solutions must be typed within this R Markdown document. When you have completed the assignment,

  1. add your name above;
  2. knit to an html file;
  3. save (or “print”) to a .pdf document; and
  4. submit your work to Moodle for grading.

Exercise 1

To linearize the discrete logistic growth model, we must take the derivative of \(f\) with respect to \(x\) then evaluate this at the equilibrium points \(x^* = 0\) and \(x^* = K\). The derivative of \(f\) is \[f'(x) = 1 + r(K - x) - rx\] Now, we evaluate the derivative at the equilibrium points: \[\text{For } x^* = 0, f'(0) = 1 + r(K-0) - r(0) = 1 + rK\] \[\text{For } x^* = K, f'(K) = 1 + r(K - K) - rK = 1 - rK\] For an equilibrium value to be stable in a discrete dynamical system, \(|f'(x^*)| < 1\). We apply this condition to the expressions found above: \[\text{For } x^*=0, \text{ stability requires } |1+rK| < 1\] \[\text{For } x^*=K, \text{ stability requires } |1-rK| < 1\] Using these inequalities, we can determine the conditions \(r\) and \(K\) for the system to be stable at each equilibrium point. \[\text{For } x^*=0, |1+rK| < 1 \text{ can be rewritten as } -1 < 1 + rK < 1\] Since \(K > 0\), safely assuming that the carrying capacity of the population is positive, we can say that \[-1<1+rK<1 = -2 < rK < 0 = \frac{-2}{K} < r < 0\] Hence, the equilibrium point \(x^* = 0\) is stable when \(r\) is in the range \(\frac{-2}{K} < r < 0\). \[\text{For } x^*=K, |1-rK| < 1 \text{ can be rewritten as } -1 < 1 - rK < 1\] Since \(K > 0\), safely assuming that the carrying capacity of the population is positive, we can say that \[-1<1-rK<1 = -2 < -rK < 0 = 0 < rK < 2 = 0 < r < \frac{2}{K}\] Hence, the equilibrium point \(x^* = K\) is stable when \(r\) is in the range \(0 < r < \frac{2}{K}\).

Exercise 2

  1. In the original model without harvesting, the equilibrium points represent values at which the population does not change from one time step to the next, i.e. \(x_{t+1} = x_t\). Adding a constant harvesting rate \(H\) introduces a constant reduction in the population at ecah time step, altering the equilibrium. The new equilibrium points \(x^*\) can be found by solving: \[x^* = x^* + rx^*(K - x^*) - H\] \[0 = rx^*(K-x^*) - H\] \[0 = -r(x^*)^2 + rKx^* - H\] \[x^* = \frac{rK \pm \sqrt{r^2K^2 - 4rH}}{2r}, x^* = 0\] are the new equilibrium points.

We can define the criteria for overfishing as the rate of harvesting that causes the fish population to decline below a certain threshold necessary for the population to sustain itself or the rate of harvesting that does not allow the population to recover to its carrying capacity in the long term. In other words, overfishing is any harvesting rate \(H\) that causes the population to go below the initial population \(x_0\) over time or does not allow it to stabilize at or above \(x_0\), or when the population’s growth becomes negative, i.e. when \(x_{t+1} < x_t\), given the harvesting rate \(H\). For experimentation purposes, we will define the threshold of overfishing as above 50% of the carrying capacity.

# Parameters
K <- 1000
r <- 0.001
x0 <- 500
threshold <- K * 0.5
time_steps <- 100

# Simulate the fish population over time with harvesting
simulate <- function(K, r, x0, H, time_steps) {
  x <- numeric(time_steps)
  x[1] <- x0
  for (t in 2:time_steps) {
    x[t] <- x[t - 1] + r * x[t - 1] * (K - x[t - 1]) - H
    if (x[t] < 0) x[t] <- 0  # Ensure population doesn't go negative
  }
  return(x)
}

# Find the threshold value of H where overfishing occurs
find_H_threshold <- function(K, r, x0, threshold, time_steps) {
  H <- 0
  overfishing_occurred <- FALSE
  while (!overfishing_occurred) {
    population <- simulate(K, r, x0, H, time_steps)
    if (any(population < threshold)) {
      overfishing_occurred <- TRUE
    } else {
      H <- H + 1
    }
  }
  return(H)
}

# Determine the threshold H
H_threshold <- find_H_threshold(K, r, x0, threshold, time_steps)

# Print the threshold value of H
print(paste("Threshold value of H for overfishing:", H_threshold))
## [1] "Threshold value of H for overfishing: 251"
# Parameters
K <- 1000
r <- 0.001
x0 <- 500
H <- 10  # Reasonable harvesting rate (1% of K)
time_steps <- 35

# Simulate the fish population over time with the chosen harvesting rate
x1 <- numeric(time_steps)
x1[1] <- x0
for (t in 2:time_steps) {
  x1[t] <- x1[t-1] + r * x1[t-1] * (K - x1[t-1]) - H
  if (x1[t] < 0) x1[t] <- 0  # Ensure the population doesn't go negative
}

# Output the result for the first 35 time steps
x1
##  [1] 500.0000 740.0000 922.4000 983.9782 989.7433 989.8948 989.8979 989.8979
##  [9] 989.8979 989.8979 989.8979 989.8979 989.8979 989.8979 989.8979 989.8979
## [17] 989.8979 989.8979 989.8979 989.8979 989.8979 989.8979 989.8979 989.8979
## [25] 989.8979 989.8979 989.8979 989.8979 989.8979 989.8979 989.8979 989.8979
## [33] 989.8979 989.8979 989.8979
# Parameters
K <- 1000
r <- 0.001
x0 <- 500
H_overharvesting <- 100  # Overharvesting rate (10% of K)
time_steps <- 35

# Simulate the fish population over time with the overharvesting rate
x2 <- numeric(time_steps)
x2[1] <- x0
for (t in 2:time_steps) {
  x2[t] <- x2[t-1] + r * x2[t-1] * (K - x2[t-1]) - H_overharvesting
  if (x2[t] < 0) x2[t] <- 0  # Ensure the population doesn't go negative
}

# Output the result for the first 35 time steps
x2
##  [1] 500.0000 650.0000 777.5000 850.4937 877.6479 885.0300 886.7819 887.1817
##  [9] 887.2720 887.2924 887.2970 887.2980 887.2983 887.2983 887.2983 887.2983
## [17] 887.2983 887.2983 887.2983 887.2983 887.2983 887.2983 887.2983 887.2983
## [25] 887.2983 887.2983 887.2983 887.2983 887.2983 887.2983 887.2983 887.2983
## [33] 887.2983 887.2983 887.2983
time_steps <- 35
time <- 1:time_steps  # Create a time vector for plotting

# Use the data from x1 and x2 to create a data frame for plotting with ggplot2
library(ggplot2)

# Create a data frame for the plot
population_data <- data.frame(
  Time = c(time, time),
  Population = c(x1, x2),
  Scenario = factor(rep(c("Reasonable Harvesting", "Overharvesting"), each = time_steps))
)

# Plot the data
ggplot(population_data, aes(x = Time, y = Population, color = Scenario)) +
  geom_line() +
  labs(title = "Fish Population Under Different Harvesting Rates",
       x = "Time (Generations)",
       y = "Population Size",
       color = "Scenario") +
  theme_minimal() +
  scale_color_brewer(palette = "Set1")

Exercise 3

# Provided data
years <- 0:34
H <- c(5577, 7450, 8130, 10307, 11012, 10850, 9519, 11183, 10521, 11767, 11351, 12749, 12465, 13132, 12028, 13679, 12349, 13193, 14136, 12700, 13378, 13731, 14043, 13719, 13399, 12131, 14769, 15327, 15489, 15279, 14626, 13152, 12945, 13698, 14293)
E <- c(279, 305, 327, 480, 589, 677, 625, 615, 586, 517, 553, 581, 648, 701, 697, 708, 785, 898, 896, 978, 1003, 964, 1043, 1163, 1096, 1099, 1168, 1333, 1851, 1905, 1858, 2307, 2303, 2334, 2305)

# Fit Model 1 using non-linear least squares
model1 <- nls(E ~ exp(a * years + b), start = list(a = 0.1, b = 5))

# Fit Model 2 using linear regression
model2 <- lm(E ~ I(years^2) + years)

# Predict values from both models
E_pred_model1 <- predict(model1, list(years = years))
E_pred_model2 <- predict(model2, list(years = years))

# Plot the actual data and the two models
library(ggplot2)
data <- data.frame(years, E, E_pred_model1, E_pred_model2)
gg <- ggplot(data, aes(years)) +
  geom_point(aes(y = E), colour = "black") +
  geom_line(aes(y = E_pred_model1, colour = "Model 1")) +
  geom_line(aes(y = E_pred_model2, colour = "Model 2")) +
  labs(title = "Fishing Effort (E) Over Time",
       x = "Time (Years)",
       y = "Fishing Effort (E)") +
  scale_colour_manual(values = c("Model 1" = "blue", "Model 2" = "red")) +
  theme_minimal() +
  theme(legend.title = element_blank())
print(gg)

  1. To determine which model is better, ww can plot the residuals of each model. The residuals of the better model should be randomly distributed and now show any pattern when plotted over time.
# Calculate residuals for Model 1
residuals_model1 <- residuals(model1)

# Calculate residuals for Model 2
residuals_model2 <- residuals(model2)

# Plot residuals for Model 1
plot(years, residuals_model1, xlab = "Years", ylab = "Residuals", main = "Residuals for Model 1")
abline(h = 0, col = "red")

# Plot residuals for Model 2
plot(years, residuals_model2, xlab = "Years", ylab = "Residuals", main = "Residuals for Model 2")
abline(h = 0, col = "red")

Based on the residuals plotted for each model, I have determined that model 2 is the better model because it displays a slightly more random distribution of the data around the horizontal line at zero on the y-axis.

effort <- function(t) {
  a <- coef(model2)["I(years^2)"]
  b <- coef(model2)["years"]
  c <- coef(model2)["(Intercept)"]
  return(a * t^2 + b * t + c)
}

Exercise 4

\[x_{t+1} = x_t + rx_t(K-x_t) - hE(t)x_t\] \[x_{t+1} = x_t + rx_t(L - x_t)\] \[hE(t)x_t = rx_t(x_t - L)\] \[L = x_t - \frac{hE(t)x_t}{rx_t}\] \[L = x_t - \frac{hE(t)}{r}\] \[L = K - \frac{hE(t)}{r}\]

\[H = hE\left(K - \frac{hE(t)}{r}\right)\] \[H = hEK - \frac{h^2E^2(t)}{r}\] \[\frac{H}{E} = hK - \frac{h^2E(t)}{r}\] Here, the term \(hK\) is the intercept and \(-\frac{h^2}{r}\) is the slope of the linear function of \(E\).

H <- c(5577, 7450, 8130, 10307, 11012, 10850, 9519, 11183, 10521, 11767, 11351, 12749, 12465, 13132, 12028, 13679, 12349, 13193, 14136, 12700, 13378, 13731, 14043, 13719, 13399, 12131, 14769, 15327, 15489, 15279, 14626, 13152, 12945, 13698, 14293)
E <- c(279, 305, 327, 480, 589, 677, 625, 615, 586, 517, 553, 581, 648, 701, 697, 708, 785, 898, 896, 978, 1003, 964, 1043, 1163, 1096, 1099, 1168, 1333, 1851, 1905, 1858, 2307, 2303, 2334, 2305)

# Calculate H/E
H_over_E <- H / E

# Plot H/E vs E
plot(E, H_over_E, xlab = "E (Effort)", ylab = "H/E (Harvesting Rate per Effort)", main = "H/E vs E")

# Example line 1
m1 <- -1/100
c1 <- 25
abline(a = c1, b = m1, col = "red")

\[h = \frac{c}{K}\] \[h = 0.025\] \[r = -\frac{h^2}{m} = - \frac{(0.025)^2}{-0.01}\] \[r = - \frac{0.000625}{-0.01}\] \[r = 0.0625\] ## Exercise 5

harvest <- function(E) {
  h <- 0.025
  r <- 0.0625
  K <- 100000
  
  L <- K - (h * E) / r
  H <- h * E * L
  return(H)
}
library(ggplot2)

# Create a data frame for plotting
data_frame <- data.frame(Year = years, Effort = E, Harvest = H)

# Predict the harvest based on the effort values
data_frame$PredictedHarvest <- sapply(data_frame$Effort, harvest)

# Plotting with ggplot2
ggplot(data_frame, aes(x = Effort)) +
  geom_point(aes(y = Harvest), color = 'blue', size = 2) +
  geom_line(aes(y = PredictedHarvest), color = 'red', size = 1) +
  labs(x = 'Effort (E)', y = 'Harvest (H)', title = 'Harvest Rate as a Function of Effort') +
  theme_minimal() +
  scale_color_manual("", 
                     breaks = c("Actual Data", "Model Prediction"),
                     values = c("Actual Data" = "blue", "Model Prediction" = "red")) +
  theme(legend.position = "bottom")

# Print out the plot
ggsave("harvest_vs_effort_plot.png", width = 10, height = 6, dpi = 300)

The model suggests a direct linear relationship between effort and harvest rate; as effort increases, the harvest rate increases proportionally. This is a simplistic model that would imply a constant return per unit of effort. However, the data points do not match this model. They are clustered at the bottom of the graph, indicating that the actual harvest rate does not increase as the effort increases, at least not to the extent predicted by the linear model.

# Parameters
r <- 0.01  # Intrinsic rate of increase
K <- 250000  # Carrying capacity
x0 <- 100000  # Initial population
t_end <- 50  # Time period for the simulation
h <- 0.0001  # Harvest rate, small value to avoid overflow

# Effort and Harvest arrays from the images
E <- c(279, 305, 327, 480, 589, 677, 625, 615, 586, 517, 553, 581, 648, 701, 697, 708, 785, 898, 896,
       978, 1003, 964, 1043, 1163, 1096, 1099, 1168, 1333, 1851, 1905, 1858, 2307, 2303, 2334, 2305)
H <- h * E

# Population model function
run_population_model <- function(x0, E, r, K, t_end, H) {
  population <- numeric(t_end + 1)
  population[1] <- x0
  
  for (t in 1:t_end) {
    H_t <- H[(t - 1) %% length(H) + 1]  # Get the harvest value for year t, cycling through the H data
    population_growth <- r * population[t] * (K - population[t])
    population[t + 1] <- population[t] + population_growth - H_t
    if (population[t + 1] < 0) {  # Population cannot be negative
      population[t + 1] <- 0
    }
  }
  return(population)
}

# Find threshold
threshold <- x0
found_threshold <- FALSE

while (!found_threshold && threshold > 0) {
  population <- run_population_model(threshold, E, r, K, t_end, H)
  if (any(population <= 0)) {
    found_threshold <- TRUE  # Population has collapsed
  } else {
    threshold <- threshold - 1000  # Decrease initial population size and try again
  }
}

# Plotting scenarios around the threshold
plot(run_population_model(threshold + 1000, E, r, K, t_end, H), type = 'l', col = 'green',
     ylim = c(0, 2.5e5), xlab = 'Time (years)', ylab = 'Population',
     main = 'Population over Time with Different Initial Sizes')
lines(run_population_model(threshold, E, r, K, t_end, H), col = 'red')
if (threshold - 1000 > 0) {
  lines(run_population_model(threshold - 1000, E, r, K, t_end, H), col = 'blue')
}
legend("topright", legend=c("Above Threshold", "At Threshold", "Below Threshold"),
       col=c("green", "red", "blue"), lty=1, cex=0.8)

# Output the threshold
print(paste("Threshold:", threshold))
## [1] "Threshold: 1e+05"

In the simulations, with a linear harvest function and the provided harvest rate (h), the long-term behavior of the species population can vary under different circumstances. With no harvesting,the population will grow logistically, approaching the carrying capacity asymptotically. This is because without harvesting, only natural growth and the environmental carrying capacity regulate the population size. With sustainable harvesting, if the harvesting function is set such that the harvest rate is sustainable, the population will also approach a stable equilibrium, which will be lower than the carrying capacity. The equilibrium is determined by the balance between the natural growth of the population and the rate at which individuals are removed from the population through harvesting.With overharvesting, if the harvesting rate is too high and unsustainable, the population could decline to extinction. Initially, there may be oscillations or a rapid decline, depending on how overharvesting is affecting the population. If the harvest rate exceeds the growth rate, the population cannot replenish itself, leading to a collapse. With threshold harvesting, when harvesting is at a threshold that the population can just withstand, the population size will stabilize at a level where growth matches the harvest rate. If the initial population size is above this threshold, the population will decline to this stable point. If the initial population size is below this threshold, it might increase to this point or collapse, depending on the exact numbers and model details. With variable harvesting, if the harvest rate varies over time, as it might in real-world scenarios, the population will show complex dynamics. It could oscillate around a certain mean level, or it might exhibit chaotic behavior, especially if the harvest rate is sensitive to the population size or other factors.

  1. From the simulations and analysis, we can draw several conclusions about how variations in the harvesting rate affect the overall fish population. The fish population is sensitive to changes in the harvesting rate. A high harvesting rate can reduce the population size significantly, potentially leading to a population collapse if it exceeds the natural replenishment rate. Setting a sustainable harvest level is crucial for the long-term stability of the fish population. Sustainable harvesting allows the population to maintain itself at a level where the birth rate balances the death rate plus the harvesting rate, ensuring that the population does not decline over time. There is a threshold effect where the population can withstand certain levels of harvesting but will collapse if that level is exceeded. Identifying this threshold is important for conservation and fishery management to avoid pushing the population to a point where it cannot recover.If the harvesting rate is managed well, the fish population can reach a dynamic equilibrium, fluctuating around a mean value that considers both the ecological capacity (carrying capacity) and the economic demands (harvesting).If the harvesting rate is not constant and fluctuates over time, the population may exhibit complex responses. These can include periodic fluctuations, potentially chaotic dynamics, or even a gradual decline if the average harvesting rate is too high over the long term. The initial population size, in relation to the carrying capacity and harvesting rate, influences how the population will respond to harvesting. A population starting at or near the carrying capacity may be more resilient to harvesting, while a smaller initial population may be more vulnerable to overharvesting. Fish population management must be adaptive, taking into account the current population size, the carrying capacity of the ecosystem, and the harvesting rate.