Python 中作为时间函数的均方位移

Question

我被分配了以下内容：

"...Plot your final MSD as a function of δt. Include errorbars σ = std(MSD)/√N, where std(MSD) is the standard deviation among the different runs and N is the number of runs. N.B.: This is the formula for the error of the mean for statistically independent, normally distributed time series. You should find that your curve resembles

MSD(δt) = Dδt,

i. e., it is linear as a function of time in spite of being the square distance the particle goes. That means that the particle undergoes diffusive motion where it advanced only proportional to the square root of time. Compute your diffusion constant D from the curve, and check that D = ∆."

但是，由于我不是很精通编码，而且我不确定我是否完全理解所涉及的计算，因此我正在努力解决这个问题。

这是我正在处理的代码，请注意，最终情节尚未完成，因为这是我正在努力处理的部分：

#%% Brownian Motion

import math
import matplotlib.pyplot as plt
import numpy as np
import random


P_dis = []

N = 10
# N = the number of iterations (i.e. the number of iterations of the loop below)
for j in range(N):
    T = 10000
    # T = time steps i.e. the number of jumps the particle makes
    x_pos = [0]
    y_pos = [0]
    # initally empty lists that will be appended below
    for i in range(1,T):
        A = random.uniform(0 , 1)
        B = A * 2 * math.pi
        current_x = x_pos[i-1]
        current_y = y_pos[i-1]
        next_x = current_x + math.cos(B)
        # takes previous xpos and applies 
        next_y = current_y + math.sin(B)
        x_pos = x_pos + [next_x]
        # appends the next x list with xpos and stores 
        y_pos = y_pos + [next_y]
    dis = np.sqrt((0 - x_pos[T - 1])**2 + (0 - y_pos[T - 1])**2)
    P_dis.append(dis)
    print(j)



plt.figure()
plt.plot(x_pos , y_pos , "0.65")

plt.plot((x_pos[0] , x_pos[-1]) , (y_pos[0] , y_pos[-1]) , "r" , label = ("Particle displacement =", dis))

plt.plot(x_pos[0] , y_pos[0] , 'ob' , label = "start" )
plt.plot(x_pos[-1] , y_pos[-1] , 'oc' , label = "end")    
plt.legend(loc = "upper left")
plt.xlabel("x position")
plt.ylabel("y position")
plt.title("Brownian Motion of a particle in 2 dimensions")
plt.grid(True)



#MSD
MSD = np.mean(P_dis)
print("Mean Square Distance is", MSD , "over" , N , "iterations")

plt.figure()
plt.plot(, [P_dis] , "r")

如有任何帮助，我们将不胜感激。

Answer 1

让我们定义：

T = 1000       # Number of time steps
N = 10         # Number of particles
step_size = 1  # Length of one step

我用 numpy 预先计算了大部分数据并将所有数据相加以获得随机游走的运动：

import numpy as np
import matplotlib.pyplot as plt

# Random direction for the N particles for T time_steps 
rnd_angles = np.random.random((N, T))*2*np.pi
# Initialize the positions for each particle to (0, 0)
pos = np.zeros((N, T, 2))                      

for t in range(1, T):
    # Calculate the position at time t for all N particles 
    # by adding a step in a random direction to the position at time t-1
    pos[:, t, 0] = pos[:, t-1, 0] + np.cos(rnd_angles[:, t]) * step_size
    pos[:, t, 1] = pos[:, t-1, 1] + np.sin(rnd_angles[:, t]) * step_size

# Calculate the distance to the center (0, 0) for all particles and all times
distance = np.linalg.norm(pos, axis=-1)

# Plot the trajectory of one particle
idx_particle = 7  # Choose from range(0, N)
x_pos = pos[idx_particle, : ,0]
y_pos = pos[idx_particle, : ,1]
dis = distance[idx_particle, -1]  # Get the distance at the last time step

plt.figure()
plt.plot(x_pos , y_pos , "0.65")
plt.plot((x_pos[0] , x_pos[-1]) , (y_pos[0] , y_pos[-1]) , "r" , 
          label=("Particle displacement =", dis))
plt.plot(x_pos[0] , y_pos[0] , 'ob' , label = "start" )
plt.plot(x_pos[-1] , y_pos[-1] , 'oc' , label = "end")
plt.legend(loc = "upper left")
plt.xlabel("x position")
plt.ylabel("y position")
plt.title("Brownian Motion of a particle in 2 dimensions")
plt.grid(True)

您可以了解正在发生的事情以及 'slow' 扩展如何通过查看一段时间内的位置来了解：

for i in np.linspace(0, T-1, 10, dtype=int):
    plt.figure()
    plt.scatter(pos[:, i, 0] , pos[:, i, 1])

您对从起点 (0, 0) 到时间的均方距离感兴趣：

squared_distance = (distance ** 2)
msd = squared_distance.mean(axis=0)
std_msd = squared_distance.std(axis=0)
sigma = std_msd / np.sqrt(N)

plt.errorbar(x=np.arange(T), y=msd, yerr=sigma)

你有机会 T, N 和 step_size 看看对msd的影响

Answer 2

要尽可能少地更改代码来绘制 MSE 及其 std，您可以执行以下操作，

import math
import matplotlib.pyplot as plt
import numpy as np
import random


P_dis = []

N = 10
# N = the number of iterations i.e. the number of iterations of the loop 
# below
T = 100
allx = np.array([]).reshape(0,T)
ally = np.array([]).reshape(0,T)
allD = np.array([]).reshape(0,T)
for j in range(N):
    # T = time steps i.e. the number of jumps the particle makes
    x_pos = [0]
    y_pos = [0]
    # initally empty lists that will be appended below
    for i in range(1,T):
        A = random.uniform(0 , 1)
        B = A * 2 * math.pi
        current_x = x_pos[i-1]
        current_y = y_pos[i-1]
        next_x = current_x + math.cos(B)
        # takes previous xpos and applies 
        next_y = current_y + math.sin(B)
        x_pos = x_pos + [next_x]
        # appends the next x list with xpos and stores 
        y_pos = y_pos + [next_y]
    dis = np.sqrt((0 - x_pos[T - 1])**2 + (0 - y_pos[T - 1])**2)
    dis = np.sqrt((0 - np.array(x_pos))**2 + (0 - np.array(y_pos))**2)
    allD = np.vstack([allD,dis])
    allx = np.vstack([allx,np.array(x_pos)])
    ally = np.vstack([ally,np.array(y_pos)])
    P_dis.append(dis)
    print(j)



plt.figure()
plt.plot(np.mean(allx,0) , np.mean(ally,0) , "0.65")
plt.figure()
plt.plot(np.arange(0,T),np.mean(allD,0),'b')
plt.plot(np.arange(0,T),np.mean(allD,0)+np.std(allD,0)/np.sqrt(N),'r')
plt.plot(np.arange(0,T),np.mean(allD,0)-np.std(allD,0)/np.sqrt(N),'r')

Python 中作为时间函数的均方位移

Mean Square Displacement as a Function of Time in Python

python

random-walk