Runge Kutta and Milstein – system of second-order coupled differential equations with noise

Question

I would like to solve a system of second-order differential equations to describe the dynamics of a system of particles. Two Newton-like forces are responsible for the motion of each particle $i$: A force acting on each particle due to other particles $f_{i,j}$ and a stochastic term of noise $f_{\mathrm{noise}}$.

The force acting on each particle due to other particles $f_{i,j}$ depends on the current position $s_i$ and velocity $v_i$ of particle $i$ and the position $s_j$ and velocity $v_j$ of the other particles $j$ of the system.

$$F_i= f_{i,j}(v_i,v_j,s_i,s_j) + f_{\mathrm{noise}}$$

The two components of motion in 2D are included for every term before mentioned. Under an Euler scheme, the velocity and position of each particle would be updated as follows: $$ \begin{alignat}{1} v_i &← v_{i} + \frac{F_i}{m} \Delta t \\ s_i &← s_i + v_i \Delta t \end{alignat} $$

where $m$ is the mass of particle $i$ and $\Delta t$ is the integration step. However, I would like to use the Milstein’s algorithm for the velocity update (since we have a term of noise) and the fourth-order Runge-Kutta method to update the position $s_i$. I am confused due to having $f_{i,j}$ dependent on $s_i$, $s_j$, $v_i$ and $v_j$. How should I operate?

Answer 1

You cannot take deterministic methods and slap a kind of Milstein correction on them to get higher order. Instead, you should be looking at methods which are specifically designed for integrating SDEs at high order if you want good accuracy. If pure efficiency is what matters and the noise is additive/diagonal, then the adaptive Rossler methods mentioned by @Wrzlprmft is the way to go. For non-diagonal, the other Rossler schemes do well.

However, in many cases when talking about "accurate" handling of second order ODEs people usually discuss symplectic integrators and energy conservation. Those methods will not conserve energy. However, there is a method due to Burrage which does conserve energy.

http://epubs.siam.org/doi/pdf/10.1137/050646032

But if you want better than strong order 0.5, then as stated in the paper, "A challenge is to extend the idea in a stochastic setting to construct explicit partitioned methods that approximate the stationary correlation matrix with high-order accuracy", i.e. you're on your own. I think a promising method to get there would be to extend Rossler's work to comparable Stochastic Runge-Kutta Nystrom methods and prove symplecticity conditions, but I don't think anyone has done that yet.

Answer 2

First of all, it is impossible to intertwine a multi-step Runge–Kutta method and the Milstein–Itō methods for a multitude of reasons that go beyond the scope of this question¹. So the best you can possibly do in is:

1. Make a deterministic Runge–Kutta step, ignoring the noise term.
2. Apply noise (in a Milstein–Itō fashion).
3. Go to 1.

This has the problem that the Runge–Kutta step totally ignores that the noise affects the system during the step. If your noise acts on much longer time scales than your Runge–Kutta steps, this may be acceptable, but then you should consider the following:

• If your Runge–Kutta time steps are this small to accurately capture the deterministic forces, i.e., the deterministic forces are acting on much shorter time scales than your stochastic forces: Do you really need noise in your model? (Note that there are certainly applications where this makes sense; just think about whether yours is one of them.)

• Otherwise, i.e., if your stochastic forces are acting on the same or even smaller time scales than your deterministic forces: Do you really gain anything from making the steps this small or use a fourth-order Runge–Kutta method in the first place? Remember that one of the main advantages of multi-step Runge–Kutta methods is to reduce the systematic errors that are the main problem of the Euler method. If you have a huge stochastic term anyway, those errors may be irrelevant to you.

That being said, there now are multi-step integrators for stochastic differential equations which combine the advantages of both sides (and are generally faster), e.g., those proposed by Rößler. Recently, Rackauckas and Nie proposed an embedded method based on this. You usually do not want to implement those yourself, but there are implementations available in Julia and Python (by yours truly).

---

^{¹ if you have a concise and digestible argument for this, please answer this question of mine}