# Predator-Prey Model with Seasonal Variation
## Overview
This model simulates the dynamics between predator and prey populations using a modified Lotka-Volterra system with seasonal variation in prey growth. The model demonstrates how environmental fluctuations can affect population cycles and ecosystem stability.
## Mathematical Model
The system is governed by the following coupled differential equations:
```
dx/dt = ax - bxy + c·sin(ωt) (prey population)
dy/dt = -dy + exy (predator population)
```
Where:
- `x(t)` = prey population at time t
- `y(t)` = predator population at time t
- `t` = time
## Parameters
### Biological Parameters
- **`a = 1.0`** - **Intrinsic prey growth rate**
- Units: 1/time
- Represents the exponential growth rate of prey in the absence of predation
- Higher values lead to faster prey population growth
- **`b = 0.5`** - **Predation rate coefficient**
- Units: 1/(predator·time)
- Determines how effectively predators consume prey
- The term `bxy` represents the rate of prey consumption based on predator-prey encounters
- **`d = 0.8`** - **Predator death rate**
- Units: 1/time
- Natural mortality rate of predators in the absence of prey
- Higher values lead to faster predator population decline when prey is scarce
- **`e = 0.3`** - **Predator efficiency**
- Units: 1/(prey·time)
- Conversion efficiency from consumed prey to new predators
- The term `exy` represents predator growth based on successful hunting
### Environmental Parameters
- **`c = 0.2`** - **Seasonal variation amplitude**
- Units: 1/time
- Magnitude of seasonal fluctuations in prey growth
- Models environmental factors like seasonal food availability, weather, or breeding seasons
- **`w = 2π/12`** - **Seasonal frequency**
- Units: radians/time
- Frequency of seasonal cycles (monthly in this case)
- The period is 12 time units, representing annual cycles with monthly resolution
## Initial Conditions
### Scenario 1: `iv1 = [4, 2]`
- Initial prey population: 4 units
- Initial predator population: 2 units
### Scenario 2: `iv2 = [4.1, 2]`
- Initial prey population: 4.1 units (slightly higher)
- Initial predator population: 2 units
- Used to demonstrate sensitivity to initial conditions
## Simulation Parameters
- **Time span**: 0 to 50 time units
- **Step size**: `h = 0.01`
- **Numerical method**: 4th-order Runge-Kutta (vectorRK4)
## Key Insights
### Equilibrium Point
Without seasonal variation (`c = 0`), the system has an equilibrium at:
- Prey equilibrium: `x* = d/e = 0.8/0.3 = 2.67`
- Predator equilibrium: `y* = a/b = 1.0/0.5 = 2.0`
### Seasonal Effects
The `c·sin(ωt)` term introduces:
- **Periodic forcing** that can synchronize population cycles
- **Resonance effects** when the forcing frequency matches natural oscillation frequency
- **Enhanced stability** or **increased chaos** depending on parameter values
### Biological Interpretation
- **Spring/Summer** (`sin(ωt) > 0`): Increased prey growth due to abundant food/favorable conditions
- **Fall/Winter** (`sin(ωt) < 0`): Reduced prey growth due to scarce resources/harsh conditions
- **Predator lag**: Predator population changes follow prey population changes with a time delay
## Output Visualizations
1. **Population Dynamics**: Time series showing both species over time
2. **Phase Plane**: Trajectory in prey-predator space showing limit cycles
3. **Sensitivity Analysis**: Comparison of trajectories with different initial conditions
4. **Parameter Effects**: Demonstrates how small changes can lead to different outcomes
## Educational Applications
This model helps students understand:
- Coupled differential equation systems
- Predator-prey relationships in ecology
- Effects of environmental variability on population dynamics
- Sensitivity analysis and chaos theory
- Numerical methods for solving ODE systems
- Phase plane analysis and limit cycles
## Codes
``` MATLAB
% Predator-Prey Model Driver
% Modified Lotka-Volterra equations with seasonal variation in prey growth
% x' = ax - bxy + c*sin(wt) (prey population)
% y' = -dy + exy (predator population)
% Parameters
a = 1.0; % intrinsic prey growth rate
b = 0.5; % predation rate
d = 0.8; % predator death rate
e = 0.3; % predator efficiency
c = 0.2; % seasonal variation amplitude
w = 2*pi/12; % seasonal frequency (monthly cycles)
% Time span
range = [0, 50];
h = 0.01;
% Initial conditions [prey, predator]
iv1 = [4, 2]; % Starting populations
iv2 = [4.1, 2]; % Slightly different initial condition
% Define the ODE system (row vector output for your functions)
predprey_ode = @(t, y) [a*y(1) - b*y(1)*y(2) + c*sin(w*t), ...
-d*y(2) + e*y(1)*y(2)];
% Solve using vectorRK4 for both initial conditions
[t1, y1] = vectorRK4(predprey_ode, range, iv1, h);
[t2, y2] = vectorRK4(predprey_ode, range, iv2, h);
% Create plots
figure;
% Time series plot
subplot(2, 2, 1);
plot(t1, y1(:, 1), 'b-', 'LineWidth', 2, 'DisplayName', 'Prey');
hold on;
plot(t1, y1(:, 2), 'r-', 'LineWidth', 2, 'DisplayName', 'Predator');
title('Population Dynamics Over Time');
xlabel('Time');
ylabel('Population');
legend;
grid on;
hold off;
% Phase plane plot
subplot(2, 2, 2);
plot(y1(:, 1), y1(:, 2), 'b-', 'LineWidth', 2, 'DisplayName', 'Trajectory 1');
hold on;
plot(y2(:, 1), y2(:, 2), 'r--', 'LineWidth', 2, 'DisplayName', 'Trajectory 2');
plot(iv1(1), iv1(2), 'bo', 'MarkerSize', 8, 'MarkerFaceColor', 'b');
plot(iv2(1), iv2(2), 'ro', 'MarkerSize', 8, 'MarkerFaceColor', 'r');
title('Phase Plane (Prey vs Predator)');
xlabel('Prey Population');
ylabel('Predator Population');
legend;
grid on;
hold off;
% Comparison of prey populations with different initial conditions
subplot(2, 2, 3);
plot(t1, y1(:, 1), 'b-', 'LineWidth', 2, 'DisplayName', 'Initial: [4, 2]');
hold on;
plot(t2, y2(:, 1), 'r--', 'LineWidth', 2, 'DisplayName', 'Initial: [4.1, 2]');
title('Prey Population Sensitivity');
xlabel('Time');
ylabel('Prey Population');
legend;
grid on;
hold off;
% Comparison of predator populations with different initial conditions
subplot(2, 2, 4);
plot(t1, y1(:, 2), 'b-', 'LineWidth', 2, 'DisplayName', 'Initial: [4, 2]');
hold on;
plot(t2, y2(:, 2), 'r--', 'LineWidth', 2, 'DisplayName', 'Initial: [4.1, 2]');
title('Predator Population Sensitivity');
xlabel('Time');
ylabel('Predator Population');
legend;
grid on;
hold off;
sgtitle('Predator-Prey Model with Seasonal Variation');
% Display equilibrium point (without seasonal variation)
eq_prey = d/e;
eq_pred = a/b;
fprintf('Equilibrium point (without seasonal variation): Prey = %.2f, Predator = %.2f\n', eq_prey, eq_pred);
```
``` Python
"""
Predator-Prey Model with Seasonal Variation
Python implementation using scipy.integrate.solve_ivp and custom RK4/Euler methods
"""
import numpy as np
import matplotlib.pyplot as plt
from scipy.integrate import solve_ivp
def euler_sys(f, t_span, y0, h):
"""
Euler's method for systems of ODEs
Parameters:
f: function that takes (t, y) and returns dy/dt
t_span: [t_start, t_end]
y0: initial conditions array
h: step size
Returns:
t: time array
y: solution array (each column is one variable)
"""
t_start, t_end = t_span
t = np.arange(t_start, t_end + h, h)
n_steps = len(t)
n_vars = len(y0)
y = np.zeros((n_steps, n_vars))
y[0] = y0
for i in range(n_steps - 1):
y[i+1] = y[i] + h * f(t[i], y[i])
return t, y
def vector_rk4(f, t_span, y0, h):
"""
4th-order Runge-Kutta method for systems of ODEs
Parameters:
f: function that takes (t, y) and returns dy/dt
t_span: [t_start, t_end]
y0: initial conditions array
h: step size
Returns:
t: time array
y: solution array (each column is one variable)
"""
t_start, t_end = t_span
t = np.arange(t_start, t_end + h, h)
n_steps = len(t)
n_vars = len(y0)
y = np.zeros((n_steps, n_vars))
y[0] = y0
for i in range(n_steps - 1):
k1 = h * f(t[i], y[i])
k2 = h * f(t[i] + h/2, y[i] + k1/2)
k3 = h * f(t[i] + h/2, y[i] + k2/2)
k4 = h * f(t[i] + h, y[i] + k3)
y[i+1] = y[i] + (k1 + 2*k2 + 2*k3 + k4) / 6
return t, y
def main():
# Parameters
a = 1.0 # intrinsic prey growth rate
b = 0.5 # predation rate
d = 0.8 # predator death rate
e = 0.3 # predator efficiency
c = 0.2 # seasonal variation amplitude
w = 2*np.pi/12 # seasonal frequency (monthly cycles)
# Time span
t_span = [0, 50]
h = 0.01
# Initial conditions [prey, predator]
iv1 = np.array([4.0, 2.0]) # Starting populations
iv2 = np.array([4.1, 2.0]) # Slightly different initial condition
# Define the ODE system
def predprey_ode(t, y):
"""
Predator-prey ODEs with seasonal variation
y[0] = prey population
y[1] = predator population
"""
dydt = np.array([
a*y[0] - b*y[0]*y[1] + c*np.sin(w*t), # prey equation
-d*y[1] + e*y[0]*y[1] # predator equation
])
return dydt
# Solve using custom RK4 for both initial conditions
t1, y1 = vector_rk4(predprey_ode, t_span, iv1, h)
t2, y2 = vector_rk4(predprey_ode, t_span, iv2, h)
# Also solve using scipy for comparison (optional)
sol1 = solve_ivp(predprey_ode, t_span, iv1, dense_output=True)
t_scipy = np.linspace(t_span[0], t_span[1], 1000)
y_scipy = sol1.sol(t_scipy).T
# Create plots
fig, axes = plt.subplots(2, 2, figsize=(12, 10))
# Time series plot
axes[0, 0].plot(t1, y1[:, 0], 'b-', linewidth=2, label='Prey')
axes[0, 0].plot(t1, y1[:, 1], 'r-', linewidth=2, label='Predator')
axes[0, 0].set_title('Population Dynamics Over Time')
axes[0, 0].set_xlabel('Time')
axes[0, 0].set_ylabel('Population')
axes[0, 0].legend()
axes[0, 0].grid(True)
# Phase plane plot
axes[0, 1].plot(y1[:, 0], y1[:, 1], 'b-', linewidth=2, label='Trajectory 1')
axes[0, 1].plot(y2[:, 0], y2[:, 1], 'r--', linewidth=2, label='Trajectory 2')
axes[0, 1].plot(iv1[0], iv1[1], 'bo', markersize=8, markerfacecolor='b')
axes[0, 1].plot(iv2[0], iv2[1], 'ro', markersize=8, markerfacecolor='r')
axes[0, 1].set_title('Phase Plane (Prey vs Predator)')
axes[0, 1].set_xlabel('Prey Population')
axes[0, 1].set_ylabel('Predator Population')
axes[0, 1].legend()
axes[0, 1].grid(True)
# Comparison of prey populations with different initial conditions
axes[1, 0].plot(t1, y1[:, 0], 'b-', linewidth=2, label='Initial: [4, 2]')
axes[1, 0].plot(t2, y2[:, 0], 'r--', linewidth=2, label='Initial: [4.1, 2]')
axes[1, 0].set_title('Prey Population Sensitivity')
axes[1, 0].set_xlabel('Time')
axes[1, 0].set_ylabel('Prey Population')
axes[1, 0].legend()
axes[1, 0].grid(True)
# Comparison of predator populations with different initial conditions
axes[1, 1].plot(t1, y1[:, 1], 'b-', linewidth=2, label='Initial: [4, 2]')
axes[1, 1].plot(t2, y2[:, 1], 'r--', linewidth=2, label='Initial: [4.1, 2]')
axes[1, 1].set_title('Predator Population Sensitivity')
axes[1, 1].set_xlabel('Time')
axes[1, 1].set_ylabel('Predator Population')
axes[1, 1].legend()
axes[1, 1].grid(True)
plt.suptitle('Predator-Prey Model with Seasonal Variation')
plt.tight_layout()
plt.show()
# Display equilibrium point (without seasonal variation)
eq_prey = d / e
eq_pred = a / b
print(f'Equilibrium point (without seasonal variation): Prey = {eq_prey:.2f}, Predator = {eq_pred:.2f}')
# Optional: Compare custom RK4 with scipy
if len(t_scipy) > 0:
plt.figure(figsize=(10, 6))
plt.plot(t1, y1[:, 0], 'b-', linewidth=2, label='Custom RK4 - Prey')
plt.plot(t1, y1[:, 1], 'r-', linewidth=2, label='Custom RK4 - Predator')
plt.plot(t_scipy, y_scipy[:, 0], 'b--', alpha=0.7, label='SciPy - Prey')
plt.plot(t_scipy, y_scipy[:, 1], 'r--', alpha=0.7, label='SciPy - Predator')
plt.title('Comparison: Custom RK4 vs SciPy solve_ivp')
plt.xlabel('Time')
plt.ylabel('Population')
plt.legend()
plt.grid(True)
plt.show()
if __name__ == "__main__":
main()
```