R中三阶ODE数值解与解析解图像不匹配问题排查

给定三阶常微分方程：$y'''(t)-3y'(t)-2y(t)=-12e^{-t}$，初始条件为$y(0)=0$，$y'(0)=1$，$y''(0)=8$。手动推导了解析解并在R中绘制，同时编写了数值解代码，但两者图像不匹配，排查出以下错误：

1. 数值解的ODE系统转换符号错误 #

原方程移项可得：
$$y'''(t) = 3y'(t) + 2y(t) - 12e^{-t}$$
但你的数值解代码中，dy3dt的表达式错误地使用了加号代替减号：

# 错误代码
dy3dt <- 3 * y2 + 2 * y1 + 12 * exp(-t)
# 正确代码
dy3dt <- 3 * y2 + 2 * y1 - 12 * exp(-t)

这是导致数值解偏离解析解的核心原因。

2. 绘图叠加逻辑问题 #

解析解的绘图代码使用了add=T，但如果数值解的绘图是单独创建的新画布（默认plot()会新建画布），add=T无法将解析解叠加到数值解的图上，导致两个曲线在不同窗口显示，看起来不匹配。需要先绘制数值解，再调用curve()叠加解析解。

3. 补充：deSolve包的加载问题 #

注意ode()函数来自deSolve包，代码中未体现加载步骤，运行前需要先安装并加载该包：

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

修正后的完整代码 #

# 加载deSolve包
install.packages("deSolve")
library(deSolve)

# 定义ODE系统（修正符号错误）
ode_system <- function(t, y, parms) {
 y1 <- y[1]  # y(t)
 y2 <- y[2]  # y'(t)
 y3 <- y[3]  # y''(t)

dy1dt <- y2
dy2dt <- y3
dy3dt <- 3 * y2 + 2 * y1 - 12 * exp(-t)

 return(list(c(dy1dt, dy2dt, dy3dt)))
}

# 初始条件
initial_conditions <- c(y1 = 0, y2 = 1, y3 = 8)

# 时间点
times <- seq(0, 5, by = 0.001)

# 求解ODE
solution <- ode(y = initial_conditions, times = times, func = ode_system, parms = NULL)
result_df <- as.data.frame(solution)

# 绘制对比图
plot(result_df$time, result_df$y1, type = 'l', col = 'blue', xlab = 'Time', ylab = 'y(t)',
 main = '解析解与数值解对比')
# 叠加解析解
O3 <- function(t) (2*t^2*exp(-t)+(2/3)*exp(2*t)-(2/3)*exp(-t)-t*exp(-t))
curve(O3, 0, 5, lty = 1, add = T, col = "green")

# 添加图例
legend("topright", legend = c("数值解", "解析解"), col = c("blue", "green"), lty = 1)

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