I am struggling with the following problem: I have a system of first-order ODEs describing a motion of particle in a give potential (actually it's 3d but doesn't matter):
dx = v dt,
dv = f(x) dt,
plus a small stochastic component, which perturbs the velocity of the particle according to
dv = g(x,v) dt + h(x,v) dW,
where dW is the Wiener process. That is, my drift and diffusion coefficients depend on the position and velocity, but they are always small compared to the derivatives in the deterministic part of ODE. I have many different methods for solving these equations (fixed-timestep 4th order RK, variable-timestep 8th order RK, Bulirsch-Stoer, etc.), which are excellent for the deterministic part. Then I thought that I could add the small perturbation to the velocity after each timestep of deterministic ODE integrator, by computing g and h at the end of timestep and slightly modifying the velocity as
dv += g dt + sqrt(h dt) * xi,
where xi is a random variable drawn from the standard normal distribution.
This works almost fine, but there is some bias in the long-term behaviour of particle energies (I compare it with the grid-based solution of Fokker-Planck equation describing the same process deterministically using the distribution function). What's worse, this bias is different (even changes sign) for different choices of ODE integrator, and doesn't go away when I decrease the timestep (actually becomes worse).
Reading various literature (e.g. Kloeden&Platen's book, or Higham(2001) introductory paper), it seems that things are much more subtle than I imagined, and I'm not even sure whether this task can be done in the way I attempt (that is, interleave the deterministic and stochastic part as in a symplectic integration kick-drift scheme). Any suggestions?