如何用Python求解含Ornstein-Uhlenbeck过程的随机Lorenz微分方程组？

问题背景 #

给定嵌入随机环境的Lorenz微分方程组：

x'(t) = (10+σ(t))(y(t)-x(t))
y'(t) = ρx(t)-y(t)-x(t)z(t)
z'(t) = x(t)y(t)-βz(t)

其中ρ、β为常数，σ(t)是Ornstein-Uhlenbeck过程，解析解为：
$$\sigma_t = \sigma_0 e^{-\alpha t} + \gamma e^{-\alpha t} \int_0^t e^{\alpha s} dW_s$$
你已通过以下Python代码完成σ(t)的样本路径模拟：

import numpy as np

t_max = 1
n = 1000
t, dt = np.linspace(0, t_max, n, endpoint=False, retstep=True)
dW = np.sqrt(dt) * np.random.randn(n)
alpha = 1
gamma = 1
sigma_0 = 0  # 需自行指定初始值
sigma = sigma_0 * np.exp(-alpha * t) + gamma * np.exp(-alpha * t) * np.cumsum(np.exp(alpha * t) * dW)

求解思路 #

由于σ(t)的样本路径已生成（每个样本下σ(t)是确定的时变函数），该系统可视为带时变参数的常微分方程组，直接使用常规数值积分方法即可求解。下面给出两种常用方法的实现：

1. 欧拉法（一阶精度） #

欧拉法是最基础的数值积分方法，迭代逻辑简单，适合快速验证。核心迭代公式为：
$$
\begin{align*}
x_{k+1} &= x_k + (10+\sigma_k)(y_k - x_k) \cdot dt \
y_{k+1} &= y_k + (\rho x_k - y_k - x_k z_k) \cdot dt \
z_{k+1} &= z_k + (x_k y_k - \beta z_k) \cdot dt
\end{align*}
$$

完整代码实现：

import numpy as np

# 设定参数与初始条件
rho = 28  # 示例值，可按需修改
beta = 8/3
sigma_0 = 0
alpha = 1
gamma = 1
t_max = 1
n = 1000
t, dt = np.linspace(0, t_max, n, endpoint=False, retstep=True)

# 生成σ(t)样本路径
dW = np.sqrt(dt) * np.random.randn(n)
sigma = sigma_0 * np.exp(-alpha * t) + gamma * np.exp(-alpha * t) * np.cumsum(np.exp(alpha * t) * dW)

# 初始化状态变量
x = np.zeros(n)
y = np.zeros(n)
z = np.zeros(n)
x[0], y[0], z[0] = 1.0, 1.0, 1.0  # 初始值可修改

# 欧拉法迭代求解
for k in range(n-1):
    dx = (10 + sigma[k]) * (y[k] - x[k]) * dt
    dy = (rho * x[k] - y[k] - x[k] * z[k]) * dt
    dz = (x[k] * y[k] - beta * z[k]) * dt
    x[k+1] = x[k] + dx
    y[k+1] = y[k] + dy
    z[k+1] = z[k] + dz

2. 龙格-库塔法（RK4，四阶精度） #

RK4是精度更高的经典数值积分方法，适合对结果精度要求较高的场景。核心迭代逻辑为计算四个中间斜率，再加权平均更新状态：

完整代码实现：

import numpy as np

# 设定参数与初始条件
rho = 28
beta = 8/3
sigma_0 = 0
alpha = 1
gamma = 1
t_max = 1
n = 1000
t, dt = np.linspace(0, t_max, n, endpoint=False, retstep=True)

# 生成σ(t)样本路径
dW = np.sqrt(dt) * np.random.randn(n)
sigma = sigma_0 * np.exp(-alpha * t) + gamma * np.exp(-alpha * t) * np.cumsum(np.exp(alpha * t) * dW)

# 初始化状态变量
x = np.zeros(n)
y = np.zeros(n)
z = np.zeros(n)
x[0], y[0], z[0] = 1.0, 1.0, 1.0

# RK4迭代求解
for k in range(n-1):
    # 计算四个中间步长
    k1_x = (10 + sigma[k]) * (y[k] - x[k]) * dt
    k1_y = (rho * x[k] - y[k] - x[k] * z[k]) * dt
    k1_z = (x[k] * y[k] - beta * z[k]) * dt

    k2_x = (10 + sigma[k]) * ((y[k] + k1_y/2) - (x[k] + k1_x/2)) * dt
    k2_y = (rho*(x[k]+k1_x/2) - (y[k]+k1_y/2) - (x[k]+k1_x/2)*(z[k]+k1_z/2)) * dt
    k2_z = ((x[k]+k1_x/2)*(y[k]+k1_y/2) - beta*(z[k]+k1_z/2)) * dt

    k3_x = (10 + sigma[k]) * ((y[k] + k2_y/2) - (x[k] + k2_x/2)) * dt
    k3_y = (rho*(x[k]+k2_x/2) - (y[k]+k2_y/2) - (x[k]+k2_x/2)*(z[k]+k2_z/2)) * dt
    k3_z = ((x[k]+k2_x/2)*(y[k]+k2_y/2) - beta*(z[k]+k2_z/2)) * dt

    k4_x = (10 + sigma[k]) * ((y[k] + k3_y) - (x[k] + k3_x)) * dt
    k4_y = (rho*(x[k]+k3_x) - (y[k]+k3_y) - (x[k]+k3_x)*(z[k]+k3_z)) * dt
    k4_z = ((x[k]+k3_x)*(y[k]+k3_y) - beta*(z[k]+k3_z)) * dt

    # 更新状态
    x[k+1] = x[k] + (k1_x + 2*k2_x + 2*k3_x + k4_x)/6
    y[k+1] = y[k] + (k1_y + 2*k2_y + 2*k3_y + k4_y)/6
    z[k+1] = z[k] + (k1_z + 2*k2_z + 2*k3_z + k4_z)/6

关键注意点 #

• σ(t)是随机过程，每次运行代码会生成不同的样本路径，对应方程组的解也会存在随机性，需多次模拟分析统计特性。
• 初始值$x(0), y(0), z(0)$以及常数ρ、β需根据具体问题需求调整。
• 若需更高精度或自适应步长，可使用scipy.integrate.solve_ivp等专业数值积分工具，只需将σ(t)作为时变参数传入右端函数即可。

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