0
$\begingroup$

I am trying to build my own simulator of Langevin Equation for the Brownian motion.

According to this material.

The way we calculate the particle position in certain time step is :

enter image description here

W(u) is a Wiener process. x0, v0, tauB are all constant.

My question is: How to write a c++ code for this Wiener process integral equation?

This is the code that I currently used in the simulator.

N=1000;
tau=0.1;      //0.1s
t=N*tau;
tauB=m/gamma; //~ 1e-8s

for(int j=0;j<=N;j++){  // t step
        sum=0;
    for(int k=0;k<=j;j++){  // u step
        sum=sum+(1-exp(-(j-k)*tau/tauB))*ND[j]*tau
    }
    x[j]=x0+v0*tauB(1-exp(-j*tau/tauB))+tauB/m*sum
}

x[j] is the particle position at time step j. ND[j] is a Normal distribution random value at time step j.

dW(t)=dU(t)=ND(t)dt.

The result is incorrect. There must be something important that I misunderstood in the equation.

Please help me.

Improve this question
$\endgroup$
  • $\begingroup$ Is your code actual code or are you using C as pseudo code? You're using square braces as parentheses in multiple places. $\endgroup$ – horchler Nov 5 '15 at 14:59
  • 1
    $\begingroup$ If your goal is to simulate this, it might be more straightforward to write out the SDE and use Euler-Maruyama. You also might find this answer of mine on Math.SE helpful: Implementing Ornstein–Uhlenbeck in Matlab. $\endgroup$ – horchler Nov 5 '15 at 15:11
  • $\begingroup$ My fault, I should write it more carefully. $\endgroup$ – Zhean Lee Nov 5 '15 at 15:12
1
$\begingroup$

The increments, $W(t_{k}) - W(t_{k-1})$, of the Brownian motion $W(t)$ are normally distributed with zero mean and variance equal to $t_{k} - t_{k-1}$. If ND[j] is a normal random variable with zero mean and variance equal to one, then you want to multiply it by sqrt(tau). The resulting ND[j] * sqrt(tau) is a normally distributed random variable with zero mean and variance equal to $t_{k} - t_{k-1}$, which is what you want.

Improve this answer
$\endgroup$
  • $\begingroup$ Ah.... yeah, that is a my typo. Let me change the description. Sorry. $\endgroup$ – Zhean Lee Nov 5 '15 at 14:49
  • 1
    $\begingroup$ I can see that you have edited the code but the symbol N coincides with the upper bound in your loop (N=1000) and, besides that, the factor tau is still present (it should not be). $\endgroup$ – Juan M. Bello-Rivas Nov 5 '15 at 14:53

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.