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.

Needs["Statistics`ContinuousDistributions`"]

Plot of the function x(t) for Δt=1

x[0] = 0

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

Δt = 1;

i = 1;

While[
      i < 11, {x[i] = 
          x[i - 1] - 
            1/ξ ksp x[
                i - 1] Δt + (Sqrt[2 Di Δt]) 
              Random[NormalDistribution[0, 1]], i = i + 1}]

plt1 = 
      ListPlot[Table[{i Δt, x[i]}, {i, 0, 10}], 
        PlotJoined ->  True, PlotStyle ->  RGBColor[1, 0, 0], 
        PlotLabel ->  "Δt=1", 
        DisplayFunction ->  Identity, Frame ->  True]

Plot of the function x(t) for Δt=0.1

x[0] = 0

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

Δt = 0.1;

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}]

plt2 = 
      ListPlot[Table[{i Δt, x[i]}, {i, 0, 100}], 
        PlotJoined ->  True, PlotStyle ->  RGBColor[1, 0, 0], 
        PlotLabel ->  "Δt=0.1", 
        DisplayFunction ->  Identity, Frame ->  True]

Plot of the function x(t) for Δt=0.01

x[0] = 0

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

Δt = 0.01;

i = 1;

While[
      i < 1001, {x[i] = 
          x[i - 1] - 
            1/ξ ksp x[
                i - 1] Δt + (Sqrt[2 Di Δt]) 
              Random[NormalDistribution[0, 1]], i = i + 1}]

plt3 = 
      ListPlot[Table[{i Δt, x[i]}, {i, 0, 1000}], 
        PlotJoined ->  True, PlotStyle ->  RGBColor[1, 0, 0], 
        PlotLabel ->  "Δt=0.01", 
        DisplayFunction ->  Identity, Frame ->  True]

Plot of the function x(t) for Δt=0.001

x[0] = 0

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

Δt = 0.001;

i = 1;

While[
      i < 10001, {x[i] = 
          x[i - 1] - 
            1/ξ ksp x[
                i - 1] Δt + (Sqrt[2 Di Δt]) 
              Random[NormalDistribution[0, 1]], i = i + 1}]

plt4 = 
      ListPlot[Table[{i Δt, x[i]}, {i, 0, 10000}], 
        PlotJoined ->  True, PlotStyle ->  RGBColor[1, 0, 0], 
        PlotLabel ->  "Δt=0.001", 
        DisplayFunction ->  Identity, Frame ->  True]

Show[GraphicsArray[{{plt1, plt2}, {plt3, plt4}}]];