Mean Square Displacement as a Function of Time in Python

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.