使用scipy.integrate将布朗运动纳入粒子轨迹积分

Question

我想在简单的线性粒子交互模型之上添加热波动.到目前为止(没有布朗运动)一切都是使用scipy.integrate.odeint完成并且工作得很好.因此,通过使用scipy.integrate metods之一找到一种包含随机运动的方法会很好.问题如下:使用Langevin热浴我必须更新粒子位置(x)和速度(v),如下所示:

x = x + v*dt

v = v +(interaction_force*dt + random_force*dt)/ mass

其中:random_force = sqrt(constant/dt)*random_number

我认为有两个问题:

1. 步骤大小dt出现在random_force中.但我不知道在运行期间通过自适应步长控制改变的当前步长.

2. 步长控制将遇到麻烦,因为只要使用两个不同的random_number来比较不同的步长,就会有相对大的差异.(我不确定是否使用了两个不同的random_numbers)

我唯一的想法是使用固定步长的方法.但我没有发现任何喷气式飞机.任何更好的想法如何解决这个问题？

Answer 1

如果您可以访问模型内​​的绝对时间，则可以使用伪随机数生成器来生成给定时间的唯一随机数。

\n\n

这种与时间相关的伪随机数生成器可以独立于步长进行操作，并且每当重复访问某个时间点时都会生成相同的值。

\n\n

如果您无法获得绝对时间，您可以简单地积分 1*dt 来获取时间。

\n\n

这是这种与时间相关的随机数生成器的非常粗略的实现：

\n\n

import numpy as np\nimport matplotlib.pyplot as plt    \nimport struct\n\n# create a gaussian "random number" that is deterministic in the time parameter\ndef t_rand(time):\n    tmp = struct.pack(\'f\', time)\n    int_time = struct.unpack(\'i\', tmp)\n    np.random.seed(int_time)\n    return np.random.randn()\n\n# demonstrate the behavior of t_rand\n\nt1 = np.linspace(0, 1, 21)\nt2 = np.linspace(0, 1, 41)\n\nplt.subplot(2, 1, 1)\nplt.plot(t1, [t_rand(t) for t in t1])\nplt.plot(t2, [t_rand(t) for t in t2])\n\nplt.subplot(2, 1, 2)\nplt.hist([t_rand(t) for t in np.linspace(0, 10, 10000)], 30)\n

\n\n

为此目的操纵全局随机种子可能并不明智，但如果问题足够孤立，它就会起作用。

\n\n

\n\n

上图显示，对于具有不同时间步长的两个序列，同一时间点会生成相同的值。下图显示了随机值的直方图，这表明至少分布的形状似乎没有受到这种对随机种子的粗暴操作的影响。

\n\n

我不知道我的方法是否真的适用于积分器。问题2可能仍然是一个问题，但我认为问题1已经得到充分解决。这可能无法完全解决问题，但可以帮助您入门。更专业的解决方案可能是使用随机微分方程求解器，例如PyS\xc2\xb3DE。

\n