含时变高斯随机变量的微分方程组如何在Mathematica中编写与求解？

问题本质与工具选择 #

你要解的是带乘性高斯噪声的随机微分方程组（SDE），普通的NDSolve默认仅支持确定性微分方程求解，不能直接处理动态随机项。高版本Mathematica的NDSolve已扩展支持SDE求解，也可以搭配内置的RandomFunction+随机过程框架实现，两种方法都适用。

随机项Subscript[\[Eta], i][t]的正确实现方式 #

不要直接在方程中插入固定的RandomReal生成的随机数，这种方式生成的随机值不会随求解步长动态更新，计算结果不符合你的模型假设。你需要将Subscript[\[Eta], i][t]对应为随机过程，按以下步骤实现：

方法1：基于ItoProcess+RandomFunction（兼容性最好） #

你原模型中Subscript[\[Eta], i][t] = m_i + σ_i * ξ_i(t)，其中ξ_i(t)为独立标准高斯白噪声，对应伊藤SDE的标准形式为：
$$dW_i(t) = \left[ (m_i - \mu - \sigma^2) W_i(t) + J(1-W_i(t)) \right] dt + \sigma_i W_i(t) dB_i(t)$$
其中$B_i(t)$为独立的标准维纳过程。

实现代码如下：

(* 自定义参数配置 *)
numb = 3; (* 方程组数量 *)
μ = 0.05; (* 全局参数 *)
globalSigmaSq = 0.02; (* 全局方差参数 *)
J = 0.1; (* 全局参数 *)
etaMeanList = {0.2, 0.3, 0.1}; (* 每个η_i的均值，可单独调整 *)
etaStdList = {0.1, 0.05, 0.2}; (* 每个η_i的标准差，可单独调整 *)
solveEnd = 10; (* 求解终止时间 *)
stepSize = 0.01; (* 采样步长 *)

(* 定义伊藤随机过程 *)
sdeProcess = ItoProcess[
  (* 漂移项（确定性部分） *)
  Table[
    D[Subscript[W, i][t], t] == (etaMeanList[[i]] - μ - globalSigmaSq)*Subscript[W, i][t] + J*(1 - Subscript[W, i][t]),
    {i, 1, numb}
  ],
  (* 待求解变量 *)
  Table[Subscript[W, i][t], {i, 1, numb}],
  (* 初始条件 *)
  Table[{Subscript[W, i], 1}, {i, 1, numb}],
  (* 时间变量 *)
  t,
  (* 独立维纳过程，对应每个η_i的独立噪声 *)
  Table[WienerProcess[], {i, 1, numb}]
];

(* 生成5条样本路径，可根据需求调整采样数量 *)
pathSamples = RandomFunction[sdeProcess, {0, solveEnd, stepSize}, 5];

(* 可视化样本路径 *)
ListLinePlot[pathSamples, PlotLegends -> Table["W_" <> ToString[i], {i, 1, numb}], Frame -> True, FrameLabel -> {"t", "W_i(t)"}]

方法2：高版本Mathematica直接用NDSolve求解 #

Mathematica 12.1及以上版本支持直接在NDSolve中指定白噪声项，写法更贴近你原始的方程格式：

sol = NDSolve[
  Join[
    (* 方程组定义 *)
    Table[
      D[Subscript[W, i][t], t] == (etaMeanList[[i]] - μ - globalSigmaSq)*Subscript[W, i][t] + J*(1 - Subscript[W, i][t]) + etaStdList[[i]]*Subscript[W, i][t]*RandomProcess`NormalWhiteNoise[t],
      {i, 1, numb}
    ],
    (* 初始条件 *)
    Table[Subscript[W, i][0] == 1, {i, 1, numb}]
  ],
  Table[Subscript[W, i], {i, 1, numb}],
  {t, 0, solveEnd},
  Method -> {"StochasticRungeKutta", "NoiseType" -> "White"}
];

(* 绘制单条求解结果 *)
Plot[Evaluate[Table[Subscript[W, i][t] /. sol, {i, 1, numb}]], {t, 0, solveEnd}, PlotLegends -> Table["W_" <> ToString[i], {i, 1, numb}]]

扩展说明 #

• 如果你需要Subscript[\[Eta], i][t]为有色高斯噪声（例如Ornstein-Uhlenbeck过程），仅需将上述代码中的WienerProcess替换为对应随机过程的定义即可。
• 随机微分方程的解为样本路径，若需要获取Subscript[W, i][t]的均值、方差等统计特征，可批量生成多条样本路径后做聚合计算。

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