在 Python 中模拟随机微分方程

问题描述

我正在尝试解决布朗粒子和朗之文动力学的 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.,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,8))
plt.plot(y,label ='Position')
plt.title("Simulation of brownian motion") 
plt.show()

结果是这样的

我想找出我的错误。请帮忙

解决方法

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

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 而不是使用 Euler 方法。这些是正确的方程式：

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