jeudi 26 mars 2020

Mean Squared Displacement for Random Walk in Python

I've been trying to compute the mean squared displacement of a random walk. The random walk function is mult_RW and outputs an array of 2d array's that are random walks starting at (0,0) that have step size 1 and step in a random direction determined by a random theta. I want to calculate the MSD of this and have tried doing so below. I'm fairly sure this is wrong but I can't figure out why as I'm fairly new to coding. Firstly, I'm not actually sure how to compute the MSD in python but my definition is:

  • MSD(δt) =〈|r(t+δt)−r(t)〉2t,ensemble

where ensemble means the mean over all the random walks I've computed and the rest of it means:

  • 〈|r(t+δt)−r(t)〉2t=(1/Nt) ∑ k|r(k+δt)−r(k)|2

My code is as follows:

def msd_calculate(N,steps):

big_array = mult_RW(N,steps) # gets our array

msd = np.zeros(steps-1) # creates an array of zeros to store the MSD for each delta*t

for i in range(1,steps-1): # iterates over each step size 
    for array in big_array: # iterates over each array
        s = [] # creates a list for the
        x, y = zip(*array) # creates seperate arrays for x and y

        x_splice, y_splice = x[1::i], y[1::i] # selects the values that are 1+i*dt
        x_splice = np.array(x_splice)
        y_splice = np.array(y_splice)
        x_diff, y_diff = np.diff(x_splice), np.diff(y_splice) # finds r(k+dt) - r(k)
        r = np.sqrt(np.square(x_diff) + np.square(y_diff)) # finds interior of sum |r(k+dt)    - r(k)|**2
        r_sq = r**2
        s.append(r_sq) #appends the value to the list so we can then calculate the mean
        '''Ignore the below commented bit I'm pretty sure that's wrong''' 
        #             r = np.sqrt(np.square(x) + np.square(y)) # calculates r for each 
        #             splice = r[1::i]
        #             r_diff = np.diff(splice)
        #             diff_sq = r_diff**2
        #             s.append(diff_sq)
    s = np.array(s)
    msd[i] = np.mean(s) # calculates mean of all of them

# Plotting the MSD function

z = [j for j in range(1,steps)] 
plt.figure(figsize=(12,12))
plt.plot(z,msd,'-b',label='MSD')
plt.xlabel(r'$\deltat')
plt.ylabel('MSD')
plt.title('MSD calculation for {} iterations of the Random Walk for {} steps'.format(N, steps))
plt.show()

return msd

If anyone can tell me whether my definition of MSD is wrong or whatever else it may be that would be amazing thank you. Also, if my code for my random walk is required (I've checked and it 100% works) I can upload that too.




Aucun commentaire:

Enregistrer un commentaire