Search code examples
pythonrandom-walk

Random walks and diffusive behaviour Python


The signature of diffusive behaviour is that the average of the square of the distance of the walker from the origin after t steps is proportional to the number of steps.
Repeat the process for several walkers and step sizes (the text suggests 500 walkers and up to 100 steps), and plot a graph similar to the one in the textbook to confirm this. Confirm also that the slope of such a graph is 1 as predicted.

I implemented the following code in Python and each time I run it the gradient I get is half the wanted value, and I could not find the mistake. Furthermore, this is the desired graph, and my graph. Can anyone spot what's wrong with the code?

import numpy as np 
import matplotlib.pyplot as plt

nwalks=500
nsteps=100
x=np.zeros(nsteps)
x2avg=np.zeros(nsteps)

array=np.zeros((nwalks,nsteps))


for j in range(nwalks):
    x2=np.zeros(nsteps)  
    for i in range(1,nsteps):
        rnd=np.random.random()
        
        if rnd<0.5:
            x[i]=x[i-1]+1
        else:
            x[i]=x[i-1]-1
        
        x2[i]=x[i]**2
        x2avg[i]=(np.sum(x2))/(i+1)
        array[j,i]=x2avg[i]

y=np.mean(array, axis=0)      
i=np.arange(nsteps)
coeffs,cov=np.polyfit(i,y,1,cov=True)

grad=coeffs[0]
intercept=coeffs[1]
dgrad=np.sqrt(cov[0][0])
dintercept=np.sqrt(cov[1][1])
print(grad,dgrad)


f1=plt.figure(1)
plt.scatter(i,y,color='black',marker='.')
#generating a function of the form y=mx + c
func = np.poly1d(coeffs)
# Getting the trendline(y values)
trendline = func(i)

plt.plot(i,trendline, 'k')

Solution

  • I have the impression that your line x2avg[i]=(np.sum(x2))/(i+1) is wrong. It calculates the mean squared distance over all steps in of the ith walker.

    Just remove the x2 and x2avg arrays and just use array[j, i] = x[i]**2 in every inner loop:

    for j in range(nwalks):
        x2=np.zeros(nsteps)
        for i in range(1,nsteps):
            rnd=np.random.random()
    
            if rnd<0.5:
                x[i]=x[i-1]+1
            else:
                x[i]=x[i-1]-1
    
            array[j,i]=x[i]**2
    

    You do already calculate the mean in the very next line, which is the correct way:

    y=np.mean(array, axis=0)