1D diffusion (Einstein & Smoluchowski)

free path length = 1; move left (-1) or right (+1) with equal probability

time
step
-8 -7 -6 -5 -4 -3 -2 -1 0 +1 +2 +3 +4 +5 +6 +7 +8 n [x2] <x2>
0 1 1 0 0
1 1 1 2 2 1
2 1 2 1 4 8 2
3 1 3 3 1 8 24 3
4 1 4 6 4 1 16 64 4
5 1 5 10 10 5 1 32 160 5
6 1 6 15 20 15 6 1 64 384 6
7 1 7 21 35 35 21 7 1 128 896 7
8 1 8 28 56 70 56 28 8 1 256 2048 8

<x2>/t = const

program diff1
! 1-dim. diffusion simulation
! Einstein & Smoluchowski
! Ulrich Schmitt, 2003-01-17

implicit none
integer, parameter :: nsteps = 8 ! number of steps
integer, dimension(-nsteps:nsteps) :: nx = 0 ! relative positions, accumulators
integer :: icount1, icount2, icount_rate, ix, irun, istep, nruns
real, dimension(nsteps) :: rnd ! random numbers
real, dimension(nsteps) :: sx=0., sx2=0. ! sum of deviations, sum of squares
real :: cpu_time

print *, 'Number of runs?'
read *, nruns

call system_clock(icount1, icount_rate) ! initial CPU time

do irun = 1, nruns
  ix = 0 ! initial position
  call random_number(rnd) ! get a sequence of random numbers
  do istep = 1, nsteps
    if (rnd(istep) < 0.5) then ! random move
      ix = ix - 1 ! left
    else
      ix = ix + 1 ! right
    end if
    sx (istep) = sx(istep) + ix
    sx2(istep) = sx2(istep) + ix**2
  end do
  nx(ix) = nx(ix) + 1 ! accumulate
end do

call system_clock(icount2) ! final CPU time
cpu_time = real(icount2 - icount1)/icount_rate ! seconds
print *, 'CPU time/s', cpu_time

print '(/, 9(1x, i6, 1x, i1))', nx
print '(/, 8(1x, f8.5))', sx/nruns ! mean deviations
print '(/, 8(1x, f8.5))', sx2/nruns ! mean squared deviations
end
 Number of runs? 256000

 CPU time/s   2.694000   (PC 150 MHz)

   1045 0   7886 0  27954 0  56149 0  69777 0  56166 0  28034 0   7997 0    992

 -0.00085 -0.00288 -0.00323 -0.00045  0.00062  0.00099  0.00305  0.00233

  1.00000  2.00255  2.99662  3.99728  4.98916  5.99683  7.00541  7.99697