如何用R语言求解非线性非自治常微分方程组？

Hey there! Solving non-autonomous, nonlinear ordinary differential equation systems in R is totally feasible, and the deSolve package is the industry standard for this kind of work—it’s robust, flexible, and handles both autonomous and non-autonomous systems seamlessly. Let me break down exactly how to implement this with a concrete example.

First, grab the package from CRAN (it’s free and well-maintained):

install.packages("deSolve")
library(deSolve)

The key rule for non-autonomous ODEs: your equation function must explicitly accept the time variable t as an input, since your system’s terms or parameters depend on time.

Let’s use a modified Lotka-Volterra system as an example—here, the prey growth rate and predator death rate fluctuate with time (the non-autonomous part), and species interactions are nonlinear:

non_autonomous_ode <- function(t, state, parameters) {
  # Unpack state variables and fixed parameters
  with(as.list(c(state, parameters)), {
    # Time-dependent terms (the non-autonomous component)
    prey_growth <- 0.6 * sin(t/8) + 0.7  # Fluctuating prey growth rate
    predator_mortality <- 0.2 * cos(t/6) + 0.15  # Fluctuating predator death rate
    
    # Nonlinear ODE equations
    d_prey <- prey_growth * prey - predation_rate * prey * predator
    d_predator <- conversion_efficiency * predation_rate * prey * predator - predator_mortality * predator
    
    # Return derivatives as a list (matches order of state variables)
    list(c(d_prey, d_predator))
  })
}

Define your starting values, fixed constants, and the time range you want to solve for:

# Initial population sizes (prey and predator)
initial_state <- c(prey = 15, predator = 3)

# Fixed parameters (predation rate, conversion efficiency)
params <- c(predation_rate = 0.04, conversion_efficiency = 0.25)

# Time sequence to solve over (from t=0 to t=120, step size 0.1)
time_points <- seq(0, 120, by = 0.1)

Use deSolve’s ode() function to run the solver. The default method (lsoda) works great for most non-autonomous nonlinear systems—it’s adaptive and handles both stiff and non-stiff problems. If you need a specific solver (like RK4 for non-stiff systems), specify it with the method argument.

# Run the solver
solution <- ode(y = initial_state, times = time_points, func = non_autonomous_ode, parms = params)

# Convert output to a data frame for easier plotting/analysis
solution_df <- as.data.frame(solution)

Plot the solution to see how your state variables change over time:

# Create a plot of prey and predator populations
plot(solution_df$time, solution_df$prey, type = "l", col = "darkgreen", 
     xlab = "Time", ylab = "Population Size", 
     main = "Solution to Non-Autonomous Nonlinear ODE System")
lines(solution_df$time, solution_df$predator, col = "darkred")
legend("topright", legend = c("Prey", "Predator"), 
       col = c("darkgreen", "darkred"), lty = 1)

• Solver Methods: If your system is stiff (some variables change extremely rapidly), try method = "bdf" or method = "lsode" instead of the default. For non-stiff systems, method = "rk4" (fourth-order Runge-Kutta) is a reliable choice.
• Debugging: If you get errors, double-check that your ODE function returns a list of derivatives in the same order as your state variables, and that all time-dependent terms correctly use the t input.
• Other Packages: While deSolve is the most widely used, the pracma package also has ODE-solving functions (like ode45), but deSolve offers more flexibility for complex systems.

内容的提问来源于stack exchange，提问作者Zen'z