I am studying the probability density function of a random walk in a confined geometry (2D-BOX). I am also comparing this probability density function to its equivalent in infinite two-dimensional plane.


The PDF for a random walk in a TWO-DIMENSIONAL BOX is obtained by multiplying the PDF for a random walk in a ONE-DIMENSIONAL BOX for the x- and y-directions. The PDF for a random walk in a ONE-DIMENSIONAL BOX is:


where $P(x,t|x_0,0)$ represents the probability that a particle move from $x_0$ to $x$ in time $t$, and $\tau = (\frac{L}{2})^2\frac{1}{D}$ , which represents the characteristic time-scale; $L$ is the length of the box (equal for each direction); $D$ is the diffusion coefficient.


In the case of a two-dimensional random walk in an infinite plane we can write:

$$P(p, t | p_0,0) = \frac{1}{4\pi Dt} \exp[-(p-p_0)^2/4Dt]$$

where $p = (x,y)$ and $p_0 = (x_0, y_0)$.

Now, my objective is to make sense of these PDFs. I implemented these functions in R (see code) and measured the PDF for time values between 1 and 10000, with time interval 100. I set D = 0.0015, L = 11, $p_0 = (2, 5)$ and $p = (8, 5)$. Since the summation to infinity cannot be execute, I chose a value of n_lim, such that the PDF does not change, which is equivalent to sum up to infinity.

prob1Dbox<-function(invL, t, invtau, x0, x, n_lim) {
  c =  pi * (pi/4) * (t * invtau)
  res = 0
  for(n in 1:n_lim){
    res = res + (exp(-1 * (n * n) * c) * cos((n * pi * x) * invL) * cos((n * pi * x0) * invL))
  return(invL + ((2 * invL) * res))

prob2Dinfinite<-function(t, r, sigma) {
  a1 = 1/(2 * pi * sigma * t)
  b1 = exp(-(r * r)/(2.0 * sigma * t))
  return(a1 * b1)

prob2DBox<-function(t,  x0,  xt_pos,  y0,  yt_pos,  invLx,  invLy, invtau_x,  invtau_y, n_lim) {
  pbx = prob1Dbox(invLx, t, invtau_x, x0, xt_pos, n_lim)
  pby = prob1Dbox(invLy, t, invtau_y, y0, yt_pos, n_lim)

If I plot the PDF of each model with respect to time, I get:

Random walk in confined geometry (2D-BOX)

enter image description here

Random walk in infinite plane

enter image description here

My concern is related to the "tail" of each distribution. In fact, in the infinite plane, we observe that the PDF decreases over time. This is because in an infinite plane a particle can indefinitely move. In the random walk in a 2D box, we observe that the PDF reaches a plateau over time. This seems to suggest that the probability to move from $p_0$ to $p$ remains constant after a certain time point. I am not sure whether it makes sense. Do you think that the PDF should be different? May be I wrongly implemented the PDF in my code.

UPDATE (29/06/2022)

Following the method of images suggested by @lpz, I implemented the proposed formula for the method of images in one-dimension as:

  prob1DboxImages<-function(t, D, x, x0_series){
  res = 0
  for(xi in x0_series){
    res = res + ((1/sqrt(4 * pi * D * t)) * exp(-((x-xi)^2)/(4 * D * t)))

I then compared the PDF obtained from prob1DboxImages with the PDF given by prob1Dbox. Using the same parameter values as in the example above, setting Zlim = 1000, Zx = c(seq(0,Zlim,1), seq(0,Zlim,1) * (-1)) and x0_series = union((xA0 + (2 * Lx * Zx)),(-xA0 + (2 * Lx * Zx))) I do get:

enter image description here

Interestingly, I do not find differences between the two PDFs, for shorter time scales, at least for the parameter values here considered.


1 Answer 1


You are correct, the boundary conditions (von Neumann along the boundary of the box), the initial conditions (Dirac delta) and the equation of motion (diffusion equation) are separable. Therefore, the solution is separable as well. Expressed in $[0,L]^2$, after expanding out the factorization: $$ P(x,y,t|x_0,y_0) = \frac{1}{L^2}\\ +\frac{4}{L^2}\sum_{n_x,n_x\in\mathbb N^*{}^2}\cos(\pi n_xx/L) \cos(\pi n_xx_0/L)\\ \cos(\pi n_yy/L) \cos(\pi n_yy_0/L)\\ e^{-(n_x^2+n_y^2)t/\tau}\\ +\frac{2}{L^2}\sum_{n_x,n_x\in\mathbb N^*\times \{0\}\cup \{0\}\times N^*}\cos(\pi n_xx/L) \cos(\pi n_xx_0/L)\\ \cos(\pi n_yy/L) \cos(\pi n_yy_0/L)\\ e^{-(n_x^2+n_y^2)t/\tau}\\ $$ with $\tau =\frac{L^2}{D\pi^2}$, which is solution to the diffusion equation $\partial_t P=D\Delta P$ with von Neuman boundary conditions. This has the advantage of showing the decay rates of the various modes.

The difference between the two is that the box is compact while the plane is not. In the first case, the walk converges exponentially to the uniform distribution $f(x)=\frac{1}{L^2}$ hence the plateau effect. You can interpret this as ergodicity, the distribution will converge to the invariant measure. You can also view this spectrally. In this case, the spectrum is discrete. At long times, only the first constant term remains, and at next to leading contribution only the two degenerate modes $n=(0,1)$ and $n=(1,0)$ will contribute and give the rate of convergence $\tau$.

In the continuous case, there is no uniform distribution so it converges to $0$. Intuitively, this is because the walk is not recurrent and random walkers run off to infinity. Furthermore, the spectrum is not gapped, so the rate of convergence is slower (power law).

Your code looks correct, and log scales will help you better test out the calculations.

Hope this helps and tell me if you need more details.


Sorry, my bad. You were right about the separability, all along, I had forgotten some modes. This time, I checked it numerically so it's correct. By the way, a good sanity check is to compare different formulas. While the above formula is good at long times, at short times, the one derived from the method of images is more accurate: $$ P(x,t|x_0) = \sum_{x'\in X}\frac{1}{\sqrt{4\pi Dt}}e^{-\frac{(x-x')^2}{4Dt}} $$ with $X = (x_0+2L\mathbb Z)\cup (-x_0+2L\mathbb Z)$.

  • $\begingroup$ Thanks @lpz, this answer is useful. However, I tried to implement your equation and compute the PDF numerically as I shown in my question. Indeed, the values I obtain using your PDF are different from doing prob1Dbox * prob1Dbox. However, with your PDF I do get negative values of probability. I'm going to update my question with a figure comparing the two approaches. Are you sure it is correct? Or may be I wrongly implemented.. $\endgroup$ Jun 29, 2022 at 13:21
  • $\begingroup$ your answer is very useful. I have few doubts. In the equation derived using the method of images, $\mathbb Z$ represents a vector of distances from the reflecting boundaries of a 1D box? For instance, let's assume that there is a particle at position $x_0 = 4$ in a 1D box of length $L = 11$, shall we define $\mathbb Z$ as the distances from x0 to the boundaries, that is $\mathbb Z = {4, 7}$ in this case? $\endgroup$ Jun 30, 2022 at 10:01
  • 1
    $\begingroup$ In the formula $\mathbb Z$ is the set of integers $\{0,1,-1,2,-2,…\}$. Forgot to scale the formula. In your case, this will give $X=\{4,-4,18,-18,26,…\}$ $\endgroup$
    – LPZ
    Jun 30, 2022 at 11:32
  • $\begingroup$ Ok, I think I have understood now. Just to clear any of my doubts, the .... in {0, 1, -1, 2, -2,..} means that the series of integers continues up to infinity? That is, $\mathbb Z = {0,1,-1,2,-2,3,-3,4,-4,5,-5,...\infty, -\infty}$. If I understood correctly, I guess that after a certain integer value, additional $x'$ values will become negligible and thus we can ignore them. right? $\endgroup$ Jun 30, 2022 at 13:56
  • 1
    $\begingroup$ Yes that’s correct. Of course for any numerical applications, you’ll truncate the series. Note that to reach the same level of accuracy, the number pf necessary terms depends strongly on $t$. For $t\to0$ the method of image expansion will need less terms, but for $t\to\infty$ the mode expansion will require less terms. In this consideration the timescale to compare $t$ is $t_0=\frac{ L^2}{2\pi D}$ $\endgroup$
    – LPZ
    Jun 30, 2022 at 16:03

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge that you have read and understand our privacy policy and code of conduct.

Not the answer you're looking for? Browse other questions tagged or ask your own question.