SDE收敛阶实验验证异常：欧拉-马里亚姆与米尔斯坦方法结果不符预期问题咨询

我正在尝试通过实验验证欧拉-马里亚姆（Euler-Maruyama）和米尔斯坦（Milstein）方法针对以下随机微分方程（SDE）的强弱收敛阶，但得到的结果与预期不符，想请教问题出在哪里。

目标SDE #

我研究的是如下耦合SDE系统（其中仅X项包含随机布朗运动，Y、Z为确定性常微分方程）：
$$
\begin{cases}
\begin{align}
dX_t&=9.5(Y_t-X_t)dt+(Y_t-X_t)dW_t\
dY_t&=(28X_t-Y_t-X_tZ_t)dt\
dZ_t&=(X_tY_t-\frac{8}{3}Z_t)dt
\end{align}
\end{cases}
$$
注：原公式中Y的方程笔误写为$Y_T$，已修正为$Y_t$

数值方法离散实现 #

我推导的两种方法离散公式如下：

欧拉-马里亚姆方法 #

$$
\begin{cases}
\begin{align}
X_{n+1}&=X_n+9.5(Y_n-X_n)\Delta t+(Y_n-X_n)\Delta W_n\
Y_{n+1}&=Y_n+(28X_n-Y_n-X_nZ_n)\Delta t\
Z_{n+1}&=Z_n+(X_nY_n-\frac{8}{3}Z_n)\Delta t
\end{align}
\end{cases}
$$

米尔斯坦方法 #

$$
\begin{cases}
\begin{align}
X_{n+1}&=X_n+9.5(Y_n-X_n)\Delta t+(Y_n-X_n)\Delta W_n-0.5(Y_n-X_n)(\Delta W_n^2-\Delta t)\
Y_{n+1}&=Y_n+(28X_n-Y_n-X_nZ_n)\Delta t\
Z_{n+1}&=Z_n+(X_nY_n-\frac{8}{3}Z_n)\Delta t
\end{align}
\end{cases}
$$

实验设置 #

• 时间区间：$[0,1]$
• 参考步长$\Delta t=3.125\times10^{-4}$：用来近似“精确解”（因为该SDE无解析解）
• 样本步长：$sample_dts = 2^5 \times dt, 2^4 \times dt, ..., 2^1 \times dt$
• 样本数量：10000次蒙特卡洛模拟
• 初始条件：$x_0=1.508870, y_0=-1.531271, z_0=25.246091$

实验代码 #

我编写的Python实现代码如下（已修正缩进问题）：

import numpy as np
import matplotlib.pyplot as plt

x0 = 1.508870
y0 = -1.531271
z0 = 25.246091
t = 1
dt = 3.125e-4
n = int(t / dt)
sample_dts = np.power(2, np.arange(5, 0, -1)) * dt
number_of_samples = 10000

def euler_maruyama(a, dt, dw):
    x, y, z = a
    return np.array([
        x + 9.5 * (y - x) * dt + (y - x) * dw,
        y + (28 * x - y - x * z) * dt,
        z + (x * y - 8 / 3 * z) * dt
    ])

def milstein(a, dt, dw):
    x, y, z = a
    return np.array([
        x + 9.5 * (y - x) * dt + (y - x) * dw - 0.5 * (y - x) * (dw ** 2 - dt),
        y + (28 * x - y - x * z) * dt,
        z + (x * y - 8 / 3 * z) * dt
    ])

def simulate(initial_condition, step_function, weak=False):
    errors = np.zeros((len(sample_dts), len(initial_condition)))
    y = np.tile(initial_condition, (1, number_of_samples))
    x = np.tile(initial_condition, (len(sample_dts), 1, number_of_samples))
    
    if not weak:
        dw = np.random.normal(scale=np.sqrt(dt), size=(n, number_of_samples))
        for i in range(n):
            y = step_function(y, dt, dw[i])
            for j, sample_dt in enumerate(sample_dts):
                scale = int(sample_dt / dt)
                if i % scale == 0:
                    x[j] = step_function(x[j], sample_dt, np.sum(dw[i * scale:(i + 1) * scale], axis=0))
        errors = np.mean(np.abs(y - x), axis=2)
    else:
        for i in range(n):
            y = step_function(y, dt, np.random.normal(scale=np.sqrt(dt), size=number_of_samples))
        for i, sample_dt in enumerate(sample_dts):
            for j in range(int(1 / sample_dt)):
                x[i] = step_function(x[i], sample_dt, np.sqrt(sample_dt) * np.random.randn(number_of_samples))
        errors = np.abs(np.mean(y, axis=1) - np.mean(x, axis=2))
    
    return errors

# 计算欧拉-Maruyama的误差
euler_maruyama_strong_errors = simulate(np.array([x0, y0, z0]).reshape((3, 1)), euler_maruyama)
euler_maruyama_weak_errors = simulate(np.array([x0, y0, z0]).reshape((3, 1)), euler_maruyama, weak=True)

# 绘制欧拉-Maruyama的误差图
for i in range(3):
    ax1 = plt.subplot(2, 1, 1)
    ax1.plot(sample_dts, euler_maruyama_strong_errors[:, i])
    ax2 = plt.subplot(2, 1, 2)
    ax2.plot(sample_dts, euler_maruyama_weak_errors[:, i])
plt.tight_layout()
plt.show()

# 计算Milstein的误差
milstein_strong_errors = simulate(np.array([x0, y0, z0]).reshape((3, 1)), milstein)
milstein_weak_errors = simulate(np.array([x0, y0, z0]).reshape((3, 1)), milstein, weak=True)

# 绘制Milstein的误差图
for i in range(3):
    ax1 = plt.subplot(2, 1, 1)
    ax1.plot(sample_dts, milstein_strong_errors[:, i])
    ax2 = plt.subplot(2, 1, 2)
    ax2.plot(sample_dts, milstein_weak_errors[:, i])
plt.tight_layout()
plt.show()

def order_of_convergence(dts, errors):
    return np.linalg.lstsq(np.column_stack((np.ones(len(dts)), np.log(dts))), np.log(errors), rcond=None)[0][1]

预期收敛阶与实际结果 #

理论预期 #

根据SDE数值方法理论：

• 欧拉-马里亚姆方法：强收敛阶为$\boldsymbol{0.5}$（即$\sqrt{\Delta t}$），弱收敛阶为$\boldsymbol{1}$（即$\Delta t$）
• 米尔斯坦方法：强收敛阶和弱收敛阶均为$\boldsymbol{1}$（即$\Delta t$），且误差幅度应明显小于欧拉-马里亚姆方法

我采用的收敛阶定义为：
$$
E{|X_T-X_N|}\leq K(\Delta t)^j\
|E{X_T}-E{X_N}|\leq K(\Delta t)^j
$$

实际计算结果 #

通过order_of_convergence函数计算得到各变量的收敛阶如下：

异常现象 #

1. 强收敛阶：欧拉-马里亚姆接近0.5但略低，米尔斯坦完全没有达到预期的1阶，且两者结果几乎一致
2. 弱收敛阶：所有结果均高于预期的1阶
3. 误差幅度：米尔斯坦方法的误差与欧拉-马里亚姆几乎相同，没有体现出精度优势

疑问 #

请问是我的Milstein离散方程推导错误，还是实验设置（比如模拟流程、误差计算方式）存在问题？

备注：内容来源于stack exchange，提问作者Midess