# Is it a correct way to calculate the Mean Square Displacement as function of time?

I would like to know if I calculate the MSD correctly, here is a Jupyter snippet I am using:

import numpy as np
import matplotlib.pyplot as plt
import numpy as np
import pandas as pd
%matplotlib inline

x = np.random.randn(100)

df = pd.DataFrame(x, columns=['x'])
df['time'] = np.arange(0,len(df))
df['MSD'] = pd.Series((x - x[0])**2).rolling(window=10, min_periods=1).mean()

x_plt, = plt.plot(x, label='x(t)')
msd_plt, = plt.plot(df['MSD'], label='MSD(t)')
plt.legend(handles=[line_up, line_down])


Is there a smart way to average over the variances? The original definition of MSD has average inside, but I don't relay averaging here because I have only one particle.

According to the article on which you base your MSD computation, it's specified as a single number (an average over N items) which is defined as the mean (average) of the sum of squared-distances between positions of $N$ particles at time $t$ and given reference positions $x_n(0)$ for each particle at a reference time.

Assuming that you want to compute MSD over a 1D particle displacement of $N$ items, you would use the following Matlab/Octave code to achieve it:

N = 1024;                 % Number of particles
std1 = 2.34;              % set some variance
xpos = std1 *sqrt(3)*2*( randn(1,N) - 0.5); % N particles with positions.
xref = 0;                 % using same reference position for each particle.
MSD = (1/N)* sum ( (xpos-xref).^2 );  % compute MSD of N particles at time t

figure,plot(xpos);        % plot positions of each particle at time t
hold on;                  % plot MSD of those N particles at time t.
plot(MSD*ones(1,N));      % MSD*ones(1,N) is used to plot N points for MSD
title('MSD over particle positions at time t');


Assuming that you also want to plot MSD as a function of sampled time, you could then then use the following code to achieve it:

N = 1024;                 % Number of particles
M = 128;                  % Number of MSD computations sample time points
std1 = 2.34;              % Set some variance, this can change over time.
xpos = std1 *sqrt(3)*2*( randn(N,M) - 0.5); % N particles for each M time samples.
xref = 0;                 % using same reference position for all time and each poarticle.
MSD = (1/N)* sum ( (xpos-xref).^2 );  % now MSD is an 1 x M vector

figure,plot(MSD);        % plot particle positions and MSD

• Why does N=1024? I am not sure how you adjusted the algorithm to one particle.
– 0x90
Jul 12, 2017 at 10:48
• N is the number of particles in the computation, you can set it according to your problem. Jul 12, 2017 at 10:51
• so for single particle N=1.
– 0x90
Jul 12, 2017 at 10:54
• yes for a single particle N=1. Jul 12, 2017 at 11:00