2

I want to simulate the underlying stochastic process of diffusion on a microscopic level and compare the result with the solution of the heat equation. However, I'm not able to match the solution of the heat equation with my computed quantities, as I'm not able to figure out the correct time scale, spatial scale and diffusion coefficient.

Consider the following algorithm in pseudo code.

Position = [0,0,0,0,...] // place particles initially at origin
                         // The format is [x1,y1,x2,y2,...]

for t=1:numberTimeSteps
   for i = 1:numberParticles
      alpha = uniform random number in [ 0, 2pi];
      steplength = normal random number with mean 0 and stdv S

      xi += sin(alpha)*steplength
      yi += cos(alpha)*steplength

   end
end

Now I have a vector of positions of all my particles. To compute the density, I count all particles that are inside a domain $[-L,L]$, which I divide into cells (i.e. I compute a histogram).

histo = zeros(Nhist,Nhist)
for i = 1 :numberParticles
    idx = SomeFunc(xi) // compute the index corresponding to position x
    idy = SomeFunc(yi) // same for y
    histo[idy][idx] +=1
end
histo = histo / numberParticles

Now I have a 2D histogram which roughly looks like the following picture. diffusion

Due to symmetrie, I now only consider the cross-section through the image. I know (or I hope to know, actually I'm not sure if this is the correct formula), that the density for a dirac impulse can be computed in 2D as \begin{align} u(x,t) = \frac{1}{(4\pi D t)^{dim/2}} exp(-x^2/(4Dt)) \end{align} where $u(x,t)$ solves the 2D heat equation on an unbounded domain $u_t = Du_{xx}$.

Now, for a large number of particles, my statistical computation should approach the analytic solution. However, I don't know how to derive $D$ from my computation. Also, I'm not sure about $x$ and $t$, as you can see in the next picture (where I apparently take the wrong values). enter image description here

So here is my question:

What analytic function with which parameters will match the computed densities from my random walk approach. How does my simulation translate into the physical framework of the heat equation.

Thank you very much!

Thomas
  • 4,363

1 Answers1

2

There is a theoretical answer to this, which appears in the derivation of the heat equation from random walks or from Brownian motion. The main result is that if you have a Brownian motion $\sigma B^x_t$ which starts at $x$, then $u(t,x)=E[\sigma B^x_t]$ solves the heat equation $u_t = \frac{\sigma}{2} u_{xx}$. So your diffusion coefficient should be $\frac{\sigma}{2}$.

You would probably be interested in checking this theoretical answer. For this you would need to choose a measure of closeness and then solve an optimization problem. The classic way would be least squares. This would be a nonlinear least squares problem, since your function depends in a nonlinear fashion on $D$. It should not be hard to solve numerically, however, especially since $\frac{\sigma}{2}$ will turn out to be a very good guess.

Ian
  • 101,645
  • Thank you very much. I assume, $\sigma$=stdv. Can you verify, if the Green's function I assume is correct? – Thomas Jun 23 '15 at 12:37
  • @sonystarmap Yes, that's right. Cf. https://en.wikipedia.org/wiki/Heat_kernel In my answer, the increments of the Brownian motion have standard deviation $\sigma \Delta t$. – Ian Jun 23 '15 at 12:38
  • And I would use unit time, i.e. $T=$Number of Time steps? – Thomas Jun 23 '15 at 12:48
  • @sonystarmap In your pseudocode, $S=\sigma \Delta t$, but you have not said what $\Delta t$ is. It is free to be chosen. – Ian Jun 23 '15 at 12:51
  • In my code, I don't really invoke $\Delta t$. Timesteps correspond to loop iterations. I will update my question with the least square fitting which looks like I'm still doing something wrong – Thomas Jun 23 '15 at 12:54
  • @sonystarmap As I said, it's free to be chosen. You can make it $1$ if you want; $D$ just rescales accordingly. (Note that in the Green's function, $D$ never appears alone, only as $Dt$.) – Ian Jun 23 '15 at 13:08
  • Okay the problem in fitting the gaussian to my computed density was in the fact, that my cells over which I counted the number of particles where not of unit size, so my solution was rescaled. Fixing this also gave me your suggested diffusion coefficient of $\sigma/2)$. – Thomas Jun 23 '15 at 15:52