# Brownian dynamics simulation

is there any way to improve computational speed of the following Brownian dynamics simulation of a trapped particle in the Markovian case ? I time-discretized the equation of motion and then use a simple Table to compute one x-value after the other.

Here is my code:

(*Initial condition*)
x[0] = 0.;
(*Constants*)
kB = 1.;
T = 1.;
gm = 1.;
kappa = 1.;
(*Time step*)
dt = 10.^-2;
(*Number of steps*)
Ntot = 10.^5;
(*Compute a single trajectory*)
Timing[
list = Table[{n*dt,
x[n] = x[n - 1] - (kappa/gm)*x[n - 1]*
dt + (1./gm) Sqrt[2. gm*kB*T]*Sqrt[dt]*
RandomReal[NormalDistribution[0., 1.]]}, {n,
Range[1, Ntot]}];]
(*Plotting the trajectory*)
ListPlot[list, Joined -> True, PlotRange -> All]


Is there a faster way to do this (regarding the Mathematica synthax/programming style)? For a high number of steps Ntot this takes a very long time unfortunately. Your help is highly appreciated. Thank you!

• You might also try ItoProcess and similar. Commented Apr 24, 2017 at 11:37

Things we can improve with small changes:

1. You are recomputing a bunch of constants at each iteration, like (1./gm), Sqrt[2. gm*kB*T]. Compute them once.

2. It is much more efficient to compute all the random numbers in one go, before the iteration starts.

3. The ideomatic way to build up lists based on previous list elements in Mathematica is FoldList, so let's try to use that.

4. This one is not about speed but style: we don't want to make capitalized variable names in Mathematica, as they might conflict with built-in symbols, and also they make the use of auto-completion more cumbersome.

Here is the modified code:

(*Constants*)
kB = 1.;
t = 1.;
gm = 1.;
kappa = 1.;
(*Time step*)
dt = 10.^-2;
(*Number of steps*)
nTot = 10.^5;

factor1 = 1 - (kappa/gm)*dt;
scaledRandomNumbers = (1./gm) Sqrt[2. gm*kB*t]*Sqrt[dt] * RandomVariate[NormalDistribution[0., 1.], nTot];

AbsoluteTiming[
list = Rest @ FoldList[{#2*dt, #1[[2]]*factor1 + scaledRandomNumbers[[#2]]}&, {0, 0.}, Range[nTot]];
]


## Multiple trajectories

Only a small modification is needed:

nTraj = 5;
scaledRandomNumbers = (1./gm) Sqrt[2. gm*kB*t]*Sqrt[dt]*
RandomVariate[NormalDistribution[0., 1.], {nTot, nTraj}];
AbsoluteTiming[
list = Rest@
FoldList[{#2*dt, #1[[2]]*factor1 +
scaledRandomNumbers[[#2]]} &, {0, ConstantArray[0., nTraj]},
Range[nTot]];]


Now list[[m]] will look like

{m*dt, {x[1][m], x[2][m], ..., x[nTraj][m]}}


To plot e.g. trajectory nr. 3 we can do

ListPlot[
Transpose[{list[[All, 1]], list[[All, 2, 3]]}]
,Joined -> True, PlotRange -> All]

• First of all thank you very much for your help. I do not just yet understand the functionality of FoldList because I have never used it before. I hope the documentation will be helpful here. I recognized that if you at another "." to the first "0" in the last line it is even faster. Commented Apr 24, 2017 at 12:02
• A look into the documentation was helpful. This improves the speed already a lot thank you. If you have any further advice to improve it, don't hesitate and let me know. Thank you! Commented Apr 24, 2017 at 16:41
• I now would like to compute several trajectories lets say nTraj=10^2. I use ParallelTable to store the random numbers in scaledRandomNumbers[m]. Then I introduce another table where I define list[m] where m displays the number of trajectorie. However, these tables are slowing it down again. Is there any smarter way how to generalize the code for several trajectories? Commented Apr 25, 2017 at 14:13
• See my edit. Looks like it's still pretty fast even for nTraj = 100. Commented Apr 25, 2017 at 14:32
• Note: I think the distribution of the random numbers is still correct in the multi-trajectory case, but I'm not an expert on that. Doing Transpose@Table[RandomVariate[NormalDistribution[0., 1.], nTot], {nTraj}] should be safe on that account though. Commented Apr 25, 2017 at 14:38