在 Python 中模拟随机微分方程

Question

我正在尝试解决布朗粒子和朗格文动力学的 SDE。 起初我尝试用普通随机数生成器模拟 2D 布朗运动， 代码是：

import numpy as np
import matplotlib.pyplot as plt

%matplotlib inline
dt = .001  # Time step.
T = 2.  # Total time.
n = int(T / dt)  # Number of time steps.
t = np.linspace(0., T, n)  # Vector of times.
sqrtdt = np.sqrt(dt)
y = np.zeros(n)
x = np.zeros(n)


for i in range(n-1):
  x[i + 1] = x[i] +  np.random.normal(0.0,1.0)
  y[i + 1] = y[i] +  np.random.normal(0.0,1.0)


fig, axs = plt.subplots(1, 1, figsize=(12, 12))
plt.plot(y, x, label ='Position')
plt.title("Simulation of Brownian motion") 
plt.show()

现在，当我尝试借助正向欧拉法模拟相同的过程时，控制方程为

m^dv/_dt=η

使用以下代码，

import numpy as np
import matplotlib.pyplot as plt


%matplotlib inline
dt = .001  # Time step.
T = 2.  # Total time.
n = int(T / dt)  # Number of time steps.
t = np.linspace(0., T, n)  # Vector of times.
sqrtdt = np.sqrt(dt)
v_x = np.zeros(n)
v_y = np.zeros(n)

y = np.zeros(n)
x = np.zeros(n)
for i in range(n-1):
  v_x[i + 1] = v_x[i] +  sqrtdt * np.random.normal(0.0,1.0)
  v_y[i + 1] = v_y[i] +  sqrtdt * np.random.normal(0.0,1.0)
  x[i+1] = x[i] + (v_x[i]*dt)
  y[i+1] = y[i] + (v_y[i]*dt)

fig, axs = plt.subplots(1, 1, figsize=(12, 8))
plt.plot(y, x, label ='Position')
plt.title("Simulation of Brownian motion") 
plt.show()

结果是这样的，

我想弄清楚我的错误。请帮忙

Answer 1

嗯，这不是一个真正的编程问题。这些行

for i in range(n-1):
  v_x[i + 1] = v_x[i] +  sqrtdt * np.random.normal(0.0,1.0)
  v_y[i + 1] = v_y[i] +  sqrtdt * np.random.normal(0.0,1.0)
  x[i+1] = x[i] + (v_x[i]*dt)
  y[i+1] = y[i] + (v_y[i]*dt)

根本不正确，因为它是 SDE。

方程的一般形式是dx = a(t, x)dt + b(t, x)dW，其中a(t, x) 是确定性的，b(t, x) 本质上是随机的（维纳过程）。把它变成数字就变成了

x[n+1] = x[n] + dx = x[n] + a(t, x[n])dt + b(t, x[n]) sqrt(dt) ξ，其中 ξ 服从均值为 0，方差为 1 的正态分布。sqrt(dt) 来自维纳过程的性质。

您应该选择Euler-Maruyama，而不是使用欧拉方法。这些是正确的方程式：

for i in range(n - 1):
    x[i + 1] = x[i] + b_x(t, x) * sqrtdt * np.random.normal(0.0, 1.0)
    y[i + 1] = y[i] + b_y(t, y) * sqrtdt * np.random.normal(0.0, 1.0)

在你的情况下b_x(t, x) = b_y(t, y) = 1