Incorporating Brownian motion into particle path integration using scipy.integrate

Question

I want to add thermal fluctuations on top of a simple linear particle interaction model. So far (without Brownian motion) everything has been done using scipy.integrate.odeint and worked fine. Therefore, it would be nice to find a way to enable random movement using one of the scipy.integrate methods. The problem is this: using the Langevin bath of the thermal bath I will have to update the particle positions (x) and velocity (v) as follows:

x = x + v * dt

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

where: random_force = sqrt (constant / dt ) * random_number

I think there are two problems:

• Step size dt appears inside random_force. But I do not know the current step size, which changes at runtime using adaptive step size control.

• The step size control will get into the problem, because as long as two different random_numbers are used to compare different step sizes, there will be a relative big difference. (I'm not sure if two different random numbers are used)

The only idea I have is to use a method with a fixed step size. But I did not find a single plane. What is the best solution, how to solve this problem?

Answer 1

If you have access to absolute time inside the model, you can use a pseudo-random number generator that generates a unique random number for a given time.

Such a time-dependent pseudo-random number generator can operate regardless of the step size and generates the same value whenever a time point is visited again.

If you do not have access to absolute time, you can simply integrate 1 * dt to get the time.

Here is a very crude implementation of such a time-dependent random number generator:

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

It is wrong to manipulate the global random seed for this purpose, but if the problem is sufficiently isolated, it will work.

The image above shows that for two sequences with different time steps, the same value is generated for the same time point. The figure below shows a histogram of random values, which indicates that at least the shape of the distribution seems untouched by this rough manipulation of random seeds.

I do not know if my approach works inside the integrator. Problem 2 can still be a problem, but I think that problem 1 is adequately resolved. This may not solve the problem completely, but it may help you get started. A more professional solution would be to use a Stochastic differential equation solver such as PyS³DE .