5
$\begingroup$

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!

$\endgroup$
10
$\begingroup$

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]
$\endgroup$
  • $\begingroup$ 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. $\endgroup$ – Display Name Apr 24 '17 at 12:02
  • $\begingroup$ 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! $\endgroup$ – Display Name Apr 24 '17 at 16:41
  • $\begingroup$ 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? $\endgroup$ – Display Name Apr 25 '17 at 14:13
  • $\begingroup$ See my edit. Looks like it's still pretty fast even for nTraj = 100. $\endgroup$ – Marius Ladegård Meyer Apr 25 '17 at 14:32
  • $\begingroup$ 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. $\endgroup$ – Marius Ladegård Meyer Apr 25 '17 at 14:38

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.