如何在R中用deSolve包数值求解含x的ODE及设置初始条件

给定带边界条件的常微分方程：

y'[x] = (a - 1)y[x]/x
y[1] = ab
其中a>0，b>0，x>0，需用R语言deSolve包的ode函数数值求解，并对比解析解。

核心问题修正 #

• 处理分母中的x：模型函数的第一个参数就是当前自变量x，直接调用即可。注意原方程要求x>0，因此求解的自变量序列不能包含0，否则会触发除以0的错误。
• 初始条件设置：初始条件对应x=1处的y值，因此求解的自变量序列要从1开始，初始值init直接设为a*b即可。

完整代码实现 #

library(deSolve)

# 定义ODE模型函数
difeq <- function(x, y, parms) {
  a <- parms["a"] # 用命名索引更直观
  b <- parms["b"]
  
  # 直接使用传入的x参数计算导数
  dy <- y * (a - 1) / x
  
  return(list(dy))
}

# 设置参数与初始条件
params <- c(a = 2, b = 1)
init <- c(y = params["a"] * params["b"]) # 对应y[1] = 2*1=2
# 生成自变量序列，从1开始（符合x>0要求，匹配初始条件）
x_seq <- seq(1, 5, by = 0.1)

# 数值求解
out <- ode(y = init, times = x_seq, func = difeq, parms = params)

# 计算解析解：y(x) = a*b*x^(a-1)
analytical <- params["a"] * params["b"] * x_seq^(params["a"] - 1)

# 对比绘图
plot(out[,1], out[,2], type = "l", col = "blue", lwd = 2, 
     xlab = "x", ylab = "y(x)", main = "数值解与解析解对比")
lines(x_seq, analytical, col = "red", lty = 2, lwd = 2)
legend("topleft", legend = c("数值解", "解析解"), 
       col = c("blue", "red"), lty = c(1,2), lwd = 2)

代码说明 #

• 模型函数直接调用传入的x参数计算导数，规避了除以0的风险（自变量序列从1开始）。
• 初始条件严格匹配方程的边界条件，确保求解起点正确。
• 新增了解析解计算与对比绘图，可直观验证数值解的准确性。

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