Draft: Leapfrog Initialization

Hi, as we discussed in !944 (merged), the LeapFrog algorithm used in YADE is not properly initialized. Then, in this MR I propose a simple way of initializing the algorithm. However, although some of the tests I did show a small increase in accuracy, other tests like energy conservation show worse performance due to the elastic and kinetic energy mismatch. I will also discuss some alternatives to this problem.

Currently, YADE docs say that the algorithm used is LeapFrog, which goes like this:

\vec{v}(t + \frac{dt}{2}) = \vec{v}(t - \frac{dt}{2}) + \frac{\vec{F}(t)}{m}dt + O(dt^3)

\vec{r}(t + dt) = \vec{r}(t) + dt \vec{v}(t + \frac{dt}{2}) + O(dt^3)

However, to achieve this, the first step should be half a step backward and then one full step forward in the velocity or just one-half step forward. If this half step is not performed, the algorithm becomes Direct Euler:

\vec{v}(t + dt) = \vec{v}(t) + \frac{\vec{F}(t)}{m}dt + O(dt^2)

\vec{r}(t + dt) = \vec{r}(t) + dt \vec{v}(t + dt) + O(dt^2)

This method is less accurate than LeapFrog. But, removes the problem of having the position and the velocity at different times.

I think the options we have are to leave everything as it is, leave everything as it is but correct the docs, or implement the LeapFrog initialization. The algorithm implemented in !944 (merged) assumes we are correctly initializing, but there's a Non-LeapFrog version that is as accurate and the modification is minimal. If we were to move away from LeapFrog, we could use Velocity Verlet instead of Direct Euler. Velocity Verlet is also second order but is better than Direct Euler.

Let me know what you think about the implementation and the other options.