在R语言中以矩阵/向量形式求解常微分方程组（使用deSolve？）

你的思路非常合理，而且这种形式在处理高维方程组时会更简洁易维护。下面我就以你提到的Lotka-Volterra模型为例，一步步演示如何实现这种矩阵形式的定义与求解。

先对齐常规写法和矩阵形式的对应关系 #

首先回顾你给出的常规Lotka-Volterra方程：

dx/dt = ax + bxy
dy/dt = dxy - cy

如果把状态变量写成向量 x = [x, y]^T，我们可以把方程拆解成元素-wise乘法和矩阵乘法的组合：

• 交叉项 x*y 可以通过矩阵乘法 M %*% x 生成，再和 x 做元素-wise乘法（R中用*表示）
• 线性项 a*x 和 -c*y 可以用向量 v 和 x 做元素-wise乘法实现

对应到你的表达式 dx = x * (M%*%x) + v * x，我们只需要构造合适的矩阵M和向量v：

# 对应常规参数a=1, b=-1, c=1, d=1（经典捕食者-猎物参数）
M <- matrix(c(0, d, b, 0), nrow=2, ncol=2)  # 对角线为0，交叉项对应b和d
v <- c(a, -c)                               # 线性项系数

完整实现代码 #

我们用R中常用的deSolve包来求解方程组，步骤如下：

1. 加载依赖包 #

# 首次使用先安装
# install.packages("deSolve")
library(deSolve)

2. 定义矩阵形式的ODE函数 #

注意deSolve的ode函数要求输入格式为function(t, x, params)，所以我们把矩阵M和向量v打包到参数列表中：

lotka_volterra_matrix <- function(t, x, params) {
  M <- params$M
  v <- params$v
  # 核心计算：元素-wise乘法(*) + 矩阵乘法(%*%)
  dx <- x * (M %*% x) + v * x
  # 返回列表形式，符合ode函数要求
  return(list(dx))
}

3. 设置参数、初始条件和时间序列 #

# 常规Lotka-Volterra参数
a <- 1
b <- -1
c <- 1
d <- 1

# 构造矩阵和向量
M <- matrix(c(0, d, b, 0), nrow=2, ncol=2)
v <- c(a, -c)

# 初始状态：猎物x=2，捕食者y=1
x0 <- c(x = 2, y = 1)

# 求解时间范围：0到20，步长0.1
times <- seq(0, 20, by = 0.1)

4. 求解并可视化结果 #

# 打包参数
params <- list(M = M, v = v)

# 调用ode求解
sol <- ode(y = x0, times = times, func = lotka_volterra_matrix, parms = params)

# 查看求解结果前几行
head(sol)

# 绘制种群数量随时间变化的曲线
plot(sol, main = "Lotka-Volterra Model (Matrix Form)", lwd = 2)

验证结果一致性 #

为了确认矩阵形式的正确性，我们可以和常规写法的结果对比：

# 常规写法的ODE函数
lotka_volterra_standard <- function(t, x, params) {
  dx <- params$a*x[1] + params$b*x[1]*x[2]
  dy <- params$d*x[1]*x[2] - params$c*x[2]
  return(list(c(dx, dy)))
}

# 求解常规形式
params_standard <- list(a=a, b=b, c=c, d=d)
sol_standard <- ode(y=x0, times=times, func=lotka_volterra_standard, parms=params_standard)

# 对比结果（数值误差可忽略）
all.equal(sol[,2], sol_standard[,2])
all.equal(sol[,3], sol_standard[,3])

注意事项 #

• 务必区分R中的*（元素-wise乘法）和%*%（矩阵乘法），这是实现矩阵形式的关键
• 确保矩阵M的维度是n×n，向量v和状态变量x的长度都是n（n为方程组维度）
• 如果你的模型包含二次项（如x²），只需调整M的对角线元素即可，扩展性很强

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