A 3D Random walk simulation

Question

I was referring to the book Computational Physics by Nicholas J. Giordano and was studying the simple Random Walk model. (Pg. 169 - 173)

I consider a 3D Random Walk where the probability of moving in the positive direction is p and in the negative direction is q, in all the 3 coordinates and p + q = 1. In my model I assume:

1. $$ p_x=q_x =0.5,\space p_y=q_y=0.5 \space and \space p_z=q_z=0.5$$
2. $$l_x = l_y=l_z=\pm 1$$

where l is the step length. In my program I consider total number of walkers = N and total number of random steps taken by each walker is s. My goals are:

1. Plot a scatter plot : $$<x_s^2> \space vs. \space t$$
2. Find the slope from the regression line.

Here is my R code:

randomWalk <- function(N, s) {
 S <- 1:s
 av <- numeric(s)
 for (r in 1:N) {
 x <- 0
 y <- 0
 z <- 0
 for (i in 1:s) {
 if (runif(1)<0.5) x <- x+1
 else x <- x-1
 if (runif(1)<0.5) y <- y+1
 else y <- y-1
 if (runif(1)<0.5) z <- z+1
 else z <- z-1
 av[i] <- av[i] + (x^2+y^2+z^2)
 }
 }
 av <- av/N
 plot(S, av,xlab="No. of steps(t)", ylab="<r^2>", main="Random Walk",cex=0.4, pch=15)
 abline(lm(av~(S)-1))
 print(lm(av~(S)-1))
}
randomWalk(500, 100) 

In the code I consider N=500 and s=100.

The plot I get is enter image description here

with value of slope approximately 3. But in page 173 there is a plot with the same condition but slope approximately 1.

My question is, Is there something wrong with my code or am I misunderstanding the plot given in page 173? A review will be helpful.

• user1149: Your inner for loop is taking 3 steps at a time: one in the x direction, one in the y direction, and one in the z direction.

Answer 1

Quoting page 173:

Another obvious generalization is to allow the walker to move in three dimensions [...] For this simulation we have restricted the steps to be of unit length along either ±x, ±y, or ±z.

I think the author meant that at each time step the walker will walk in one of three random directions (x or y or z), not all three. Without making drastic changes to your code, you could replace this section:

 if (runif(1)<0.5) x <- x+1 else x <- x-1
 if (runif(1)<0.5) y <- y+1 else y <- y-1
 if (runif(1)<0.5) z <- z+1 else z <- z-1

with the following:

 direction <- sample(c("x", "y", "z"), 1)
 if (direction == "x") { if (runif(1)<0.5) x <- x+1 else x <- x-1 }
 if (direction == "y") { if (runif(1)<0.5) y <- y+1 else y <- y-1 }
 if (direction == "z") { if (runif(1)<0.5) z <- z+1 else z <- z-1 }

With this change, the slope is approximately 1, as shown in Figure 7.11.