python - 在 Python 中模拟随机微分方程

Question

我正在尝试求解布朗粒子和朗之万动力学的 SDE。起初我尝试用普通随机数生成器模拟二维布朗运动，代码是:

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) xi ，其中 xi 服从均值 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