Mean Square Displacement as a Function of Time in Python

Question

I have been assigned the following:

"...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 = ∆."

However, I am struggling with this as I am not very proficient at coding, as well as the fact that I am not sure I totally understand the calculation involved.

Here is the code I am working on, note that the final plot is unfinished as this is the part I am struggling with:

#%% 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")

Any help on the matter would be greatly appreciated.

Autolysin

Answer 1

To plot the MSE with its std with as little changes to your code as possible, you can do the following,

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')

Rice

hi, thanks for the reply, when I run this code, I get ValueError: all the input array dimensions except for the concatenation axis must match exactly – Autolysin

N was defined after using it, just changed it. Anyways, that's not the problem you report. try cleaning your environment. – Rice

Thank you so much for this, was about to put my head through my monitor – Autolysin

sorry, just to clarify, the red curves on the plot are the error, correct? – Autolysin

Yes, but I just noticed they are just the standard deviation. To get the errors you have to divide by sqrt(N). This is, in the last two lines, you have to change np.std(allD,0) for np.std(allD,0)/np.sqrt(N) – Rice

after doing that you should see that the error bars increase on time but are much smaller. – Rice

Let us continue this discussion in chat. – Autolysin

Answer 2

Let's define:

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

I precompute most of the data with numpy and add everything up to get the motion of the random walk:

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)

You can get an idea of what is happening and how 'slow' the expansion is moving on by looking at the positions over the time:

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

You are interested in the mean squared distance from the start point (0, 0) with respect to the time:

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)

You can chance T, N and step_size to look at the influence on msd.

Parliamentarian

thanks for the reply, I'm currently trying to sift through this and make sense of it. As I said, I'm really not very good at coding so I really only have a basic knowledge of python, meaning very limited in terms of what a few of these lines actually do – Autolysin

That is to say, how can I combine this with the initial part of my code which actually plots the motion of the particles? – Autolysin

Sorry, I should have stayed closer to your code, I added your plot now – Parliamentarian