含白噪声与狄拉克δ函数的ODE数值求解方法咨询

一、核心问题拆解 #

你的方程组包含两类特殊项，需针对性处理：

• 白噪声（W）：属于随机微分方程（SDE）范畴，不能用普通ODE求解器，需用SDE专用工具
• 狄拉克δ函数：本质是瞬时脉冲，可通过在指定时间点直接修改状态变量实现

二、R语言实现方案 #

1. 依赖包准备 #

需要sde包处理随机微分方程，deSolve辅助脉冲逻辑：

install.packages(c("deSolve", "sde"))
library(deSolve)
library(sde)

2. 处理狄拉克δ函数（脉冲项） #

狄拉克δ函数$\delta(t-t_0)$表示在$t=t_0$时刻给状态变量施加瞬时增量，我们用事件函数在指定时间点修改状态：

# 事件函数：在指定时间点给v施加脉冲（示例时间点可按需修改）
pulse_event <- function(t, state, parameters) {
  with(as.list(c(state, parameters)), {
    # 触发脉冲的时间点
    pulse_times <- c(10, 30, 50)
    if (t %in% pulse_times) {
      # kappa1为δ项的系数，需提前加入参数列表
      state["v"] <- state["v"] + kappa1
    }
    return(state)
  })
}

3. 处理白噪声的随机微分方程 #

将原方程组改写为伊藤型随机微分形式，白噪声对应扩散项：

# 定义SDE的漂移项（原ODE右侧）和扩散项（白噪声系数）
sde_model <- function(t, x, theta) {
  alpha <- theta[1]
  beta <- theta[2]
  kappa2 <- theta[3]
  kappa3 <- theta[4]
  
  # 漂移项：对应原ODE的dv/dt和dp/dt
  drift_v <- -2*alpha*x[1] - beta^2*x[2] + kappa3*x[1]^3
  drift_p <- x[1] + kappa2*x[1]^2
  
  # 扩散项：白噪声W的系数（sigma_v为v的噪声系数，需加入参数）
  diff_v <- theta[5]
  diff_p <- 0  # 若p的方程也有白噪声，修改此处数值
  
  return(list(c(drift_v, drift_p), matrix(c(diff_v, 0, 0, diff_p), nrow=2)))
}

4. 完整求解流程 #

# 补充参数：加入脉冲系数kappa1和白噪声系数sigma_v
parameters <- c(alpha = 1, beta = 2, kappa2 = 1, kappa3 = -1, kappa1 = 0.5, sigma_v = 0.1)
state <- c(v=0.8, p=0.5)
times <- seq(0, 100, by=0.01)

# 第一步：求解SDE得到基础随机轨迹
set.seed(123)  # 设置随机种子保证结果可复现
sde_out <- sde.sim(X0=state, t0=times[1], T=times[length(times)], delta=0.01, drift=sde_model, theta=parameters)
sde_df <- as.data.frame(sde_out)
colnames(sde_df) <- c("time", "v", "p")

# 第二步：应用脉冲事件修改轨迹
for (t in times) {
  idx <- which(sde_df$time == t)
  sde_df[idx, ] <- pulse_event(t, sde_df[idx, ], parameters)
}

# 可视化结果
par(oma = c(0, 0, 3, 0))
plot(sde_df$time, sde_df$v, type="l", xlab="time", ylab="value", col="blue")
lines(sde_df$time, sde_df$p, col="red")
legend("topright", legend=c("v", "p"), col=c("blue", "red"), lty=1)
mtext(outer = TRUE, side = 3, "SDE with Pulse", cex = 1.5)

三、关键说明 #

• 白噪声调整：若p的方程也包含白噪声，只需修改sde_model中的diff_p数值即可
• 脉冲灵活性：若脉冲时间是随机分布的，可通过抽样生成pulse_times；也可使用deSolve的events参数将脉冲逻辑整合到求解流程中
• 参数校准：需根据你的实际方程组，补充kappa1（δ项系数）和sigma_v（白噪声系数）的具体数值

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