欧拉-马尔可夫方法求解OU过程中<a[t]da[t]>数值与解析值不符问题

问题背景 #

给定奥恩斯坦-乌伦贝克过程的随机微分方程：
$$da[t] = -k a[t] dt + \sqrt{k n} w[t] dt$$
（注：式中w[t]为白噪声，满足<w(t)w(t')>=δ(t-t')）

使用欧拉-马尔可夫（Euler-Maruyama）方法求解a[t]的原代码可正确复现解析解的均值及高阶矩，但扩展计算变量与时间导数乘积的期望<a[t]·da[t]>时，数值结果（k=n=1时约为-0.4）与解析解（k=n=1时为1/2）严重不符。

错误分析 #

1. a的更新公式错误 #

扩展代码中，a的随机项系数被错误修改：

# 错误写法：多嵌套了一层sqrt(dt)
a=a+(-kappa*a)*dt+np.sqrt(dt)*np.sqrt(dt*kappa*nth)*w[n,m]

原正确代码的随机项系数为np.sqrt(dt)*np.sqrt(kappa*nth)*w[n,m]，额外的sqrt(dt)会导致随机项缩放错误，直接破坏a的分布特性。

2. da[t]的离散定义错误 #

连续SDE中的da[t]是微分，对应离散场景下的增量Δaₙ = aₙ₊₁ - aₙ，因此da[t]的合理离散近似应为Δaₙ / dt，而非直接套用SDE的右边项。原因在于：

• 连续白噪声w[t]的离散化形式是sqrt(dt)*w[n,m]（w[n,m]~N(0,1)），SDE右边的离散形式为-κaₙ + sqrt(κn)/sqrt(dt)*w[n,m]，直接用该式作为da[t]会引入量级错误，且忽略了aₙ与当前步噪声的相关性问题。

修正后的代码 #

核心修正点 #

• 恢复正确的a更新公式
• 用增量除以步长作为da[t]的离散近似
• 增加样本数、减小步长提升数值精度

import numpy as np

kappa = 1
nth = 1
Mmax = 10000  # 提升样本数量降低统计误差
Tmax = 5.0
dt = 0.05  # 减小步长优化离散近似精度
Nmax = int(Tmax / dt)

t_list = np.arange(0, Tmax + dt/2, dt)
w = np.random.randn(Nmax, Mmax)  # 仅需Nmax个噪声样本（对应Nmax个时间步）

a_list = np.zeros((Nmax+1, Mmax), dtype=np.complex64)
ex_list = np.zeros((Nmax+1, Mmax), dtype=np.complex64)

for m in range(Mmax):
    a = 0 + 1j*0
    a_list[0, m] = a
    # 从第1个时间步开始计算增量
    for n in range(Nmax):
        # 正确的Euler-Maruyama更新
        delta_a = (-kappa * a) * dt + np.sqrt(dt) * np.sqrt(kappa * nth) * w[n, m]
        a_new = a + delta_a
        # 用增量除以步长近似da[t]
        da = delta_a / dt
        # 记录当前时刻a与da的乘积
        ex_list[n, m] = a * da
        # 更新a到下一个时刻
        a = a_new
        a_list[n+1, m] = a

# 计算每个时间点的期望
ex_mean = np.array([np.mean(ex_list[n, :]) for n in range(len(t_list))])

# 打印稳态阶段的平均结果（t足够大时）
print(f"稳态<a·da>数值期望：{np.mean(ex_mean[-10:]).real:.4f}")
print(f"解析解结果：{kappa * nth / 2}")

结果验证 #

当kappa=1、nth=1时，增大样本数至10000并将步长设为0.05，数值结果会趋近于解析解0.5，解决原代码的偏差问题。

内容的提问来源于stack exchange，提问作者J.Agusti