Brownian Dynamics of a Spherical Particle
in a Quadratic Energy Well
Author's Data: Housam Binous
Department of Chemical Engineering
National Institute of Applied Sciences and Technology
Tunis, TUNISIA
Email: binoushousam@yahoo.com

Acknowledgement :

This example was given by Kenneth J. Beers in Numerical Methods for Chemical Engineering, Applications in Matlab, Cambridge University Press, 2007.

Brownian dynamics simulation

Needs["Statistics`ContinuousDistributions`"]

x[0] = 0

ξ  = 1; ksp = 1; Di = 1;

Δt = 0.01;

i = 1;

For[j = 1, 
      j < 5001, {i = 1, 
        While[i < 
            101, {x[i] = 
              x[i - 1] - 
                1/ξ ksp x[
                    i - 1] Δt + (Sqrt[2 Di Δt]) Random[NormalDistribution[0, 1]], i = i + 1}], 
        X[j] = x[i - 1], j++, 
        If[Mod[j, 100] ==  0, Print[j]]}]

Probability distribution P(x)

<< Statistics`DataManipulation`

t = BinCounts[Table[X[j], {j, 1, 4999}], {-5, 5, 0.05}]

midpoints = Table[i + 0.025, {i, -5, 5 - 0.05, 0.05}]

su = Sum[t[[i]], {i, 1, Length[t]}]*0.05

plt1 = 
      ListPlot[Table[ {midpoints[[i]], t[[i]]/su}, {i, 1, 
            Length[t]}], PlotJoined ->  True,  
        PlotStyle  ->  RGBColor[1, 0, 0]]

kbT = 1;

plt2 = 
      Plot[Exp[-x^2/(2*kbT)]/
          NIntegrate[Exp[-u^2/(2*kbT)], {u, -4, 4}], {x, -5, 
          5}, PlotStyle ->  RGBColor[0, 0, 1]]

Show[plt1, plt2]