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 :
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.