R语言Lotka-Volterra模型状态比例循环实现技术求助

Got it, let's work through this step by step. First, a quick note: your original initial state c(x = 0, y = 0) will just give you all zeros forever—you can't grow a population from nothing! We'll fix that first, then set up a loop to test different initial ratios of aphids (x) to parasitoids (y).

Step 1: Keep your core model function (it's solid!) #

First, we'll retain your Lotka-Volterra model definition since that's correct:

lotVmod <- function(Time, State, Pars) {
  with(as.list(c(State, Pars)), {
    dx <- x * (alpha - beta * y)
    dy <- y * (gamma - delta * x)
    return(list(c(dx, dy)))
  })
}

# Fixed parameters (you can tweak these later if needed)
Pars <- c(alpha = 1, beta = 1.1, gamma = 0.5, delta = 1)
Time <- seq(0, 100, by = 1)

Step 2: Define the initial ratios to test #

Let's pick a set of initial ratios (e.g., aphids making up 10% to 90% of the total initial population) and fix a total initial population size (so we're testing proportional changes, not just absolute numbers):

# Ratios of aphids (x) to total initial population
initial_ratios <- seq(0.1, 0.9, by = 0.2)
# Total initial population (adjust this to your needs)
total_initial_pop <- 10

Step 3: Loop through each ratio and run the model #

We'll use a loop to iterate over each ratio, calculate the corresponding initial x and y values, run the ODE solver, and store all results. We'll also add a label to each result set so we can track which ratio it came from:

library(deSolve)  # Make sure you have this loaded for `ode()`
library(dplyr)    # For combining results easily

# Initialize a list to store results from each ratio
results_list <- list()

for (i in seq_along(initial_ratios)) {
  current_ratio <- initial_ratios[i]
  # Calculate initial x and y based on the ratio
  State <- c(
    x = total_initial_pop * current_ratio,
    y = total_initial_pop * (1 - current_ratio)
  )
  # Run the model
  model_output <- as.data.frame(ode(func = lotVmod, y = State, parms = Pars, times = Time))
  # Add a column to track the initial ratio
  model_output$initial_ratio <- paste0("x:y = ", round(current_ratio, 1), ":", round(1 - current_ratio, 1))
  # Save to the list
  results_list[[i]] <- model_output
}

# Combine all results into a single data frame for easier plotting
combined_results <- bind_rows(results_list)

Step 4: Visualize the results #

We can use ggplot2 for clean, readable plots that let you compare all ratios at once. First, we'll reshape the data to long format:

library(ggplot2)
library(tidyr)

# Reshape data to long format (better for ggplot)
long_results <- combined_results %>%
  pivot_longer(cols = c(x, y), names_to = "Species", values_to = "Population") %>%
  mutate(Species = case_when(
    Species == "x" ~ "Aphids",
    Species == "y" ~ "Parasitoid"
  ))

# Plot all ratios and species
ggplot(long_results, aes(x = time, y = Population, color = initial_ratio, linetype = Species)) +
  geom_line(linewidth = 1) +
  labs(
    x = "Time",
    y = "Population Size",
    color = "Initial x:y Ratio",
    linetype = "Species"
  ) +
  theme_minimal() +
  theme(legend.position = "top")

If you prefer base R (like your original code), here's a version using matplot:

# Set up colors and line types for each ratio
colors <- rainbow(length(initial_ratios))
ltys <- 1:length(initial_ratios)

# Plot the first ratio's aphids to set up the plot
first_ratio <- unique(combined_results$initial_ratio)[1]
plot(
  x = combined_results$time[combined_results$initial_ratio == first_ratio],
  y = combined_results$x[combined_results$initial_ratio == first_ratio],
  type = "l",
  col = colors[1],
  lty = ltys[1],
  xlab = "Time",
  ylab = "Population",
  ylim = range(combined_results[, c("x", "y")])
)

# Add remaining aphid curves
for (i in 2:length(initial_ratios)) {
  current_ratio <- unique(combined_results$initial_ratio)[i]
  lines(
    x = combined_results$time[combined_results$initial_ratio == current_ratio],
    y = combined_results$x[combined_results$initial_ratio == current_ratio],
    col = colors[i],
    lty = ltys[i]
  )
}

# Add parasitoid curves (use thicker lines to distinguish)
for (i in 1:length(initial_ratios)) {
  current_ratio <- unique(combined_results$initial_ratio)[i]
  lines(
    x = combined_results$time[combined_results$initial_ratio == current_ratio],
    y = combined_results$y[combined_results$initial_ratio == current_ratio],
    col = colors[i],
    lty = ltys[i],
    lwd = 1.5
  )
}

# Add legend
legend(
  "topleft",
  legend = c(
    sapply(unique(combined_results$initial_ratio), function(r) paste(r, "(Aphids)")),
    sapply(unique(combined_results$initial_ratio), function(r) paste(r, "(Parasitoid)"))
  ),
  col = rep(colors, 2),
  lty = rep(ltys, 2),
  box.lwd = 0,
  cex = 0.8
)

Why your earlier loops might have failed #

• Zero initial populations: Starting with x=0 and y=0 gives no population growth, so all outputs are zero—this makes it look like the loop isn't working.
• Not storing results: If you didn't save each iteration's output to a list/data frame, you'd only see the last run's results.
• Plotting issues: Trying to plot all curves at once without distinguishing colors/line types makes them overlap unreadably.

内容的提问来源于stack exchange，提问作者Stuart