**Quote from frankuitaalst on 05/24/15 at 04:53:52:**

I mean : introducing higher order in the integration may allow to set larger timesteps for a given "accurancy" .

A simple test may be performed : set up a system with central mass and a test particle .

The particle should stay within ...x AU after T time at a timestep of dt .

Anxious to see how the integrator performs in GravSim ...

I personally use sometimes a 21-order Taylor integrator when I want accurate simulations with self adapting time step .

I tested this one-body verion with circular orbit, and there are some general results:

I just set the 2-body sys. with parameters:

y[0].m = 1; y[0].r.x = 0; y[0].v.y = 0; // the central mas, which stays only

y[1].m = 0; y[1].r.x = 1; y[1].v.y = 1; // a null particle

The parameters of the orbit are:

e = 0, a = 1, GM = 1, then a period is just:

T = 2pi (m/a^3)^0.5 = 2pi * 1 = 2pi

and a speed:

v = 2pi/T = 1

The error is measured as: er = r - 1 = (x^2+y^2)^0.5 - 1

And it is about 1e-9 - 1e-8 for 10 milions of steps, ie. it's no less for any order.

For 2-order method - it's the standard verlet actually,

the error is always rather big: 1e-4 or even more for big steps (below about 100 for one orbit: dt < 2pi/100).

For the 4-order and higher the error is 1e-10, and it grows with a time steadily.

There is some limit for a time step size for any order...

In the extremal case: order 16 it's about 10 steps for orbit... and for 20-thy even 4 is sufficient.

The symplectic method 8-order has similar error, as the 8-order extrapolation,

but the symplectic orbit is deformed: eccentricity e, and the semi-axis 'a' changes

rather strongly - up to 0.001 even;

The extrapolated 8-order version keeps perfectly the parameters of the orbit!

And the 10-order extrapolated method is high more precise already, than the symplectic 8-order.

...

The 12-order variant (k=6) looks very nice - probably the optimal one, for double precision computations.

The simplectic 8-order code is alsor very simple:

**Code:**void symplectic8()
{
const double W1 = -0.161582374150097E1;
const double W2 = -0.244699182370524E1;
const double W3 = -0.716989419708120E-2;
const double W4 = 0.244002732616735E1;
const double W5 = 0.157739928123617E0;
const double W6 = 0.182020630970714E1;
const double W7 = 0.104242620869991E1;
const double W0 = (1-2*(W1+W2+W3+W4+W5+W6+W7));
const static double d8[15] = { W7, W6, W5, W4, W3, W2, W1, W0, W1, W2, W3, W4, W5, W6, W7 };
const static double c8[16] = { W7/2, (W7+W6)/2, (W6+W5)/2, (W5+W4)/2,
(W4+W3)/2, (W3+W2)/2, (W2+W1)/2, (W1+W0)/2,
(W1+W0)/2, (W2+W1)/2, (W3+W2)/2, (W4+W3)/2,
(W5+W4)/2, (W6+W5)/2, (W7+W6)/2, W7/2 };
for(int i = 0; i < 16; i++)
{
h = dt * c8[i];
for(int j = 0; j < nB; j++)
y[j].r.addm(y[j].v, h);
if( i != 15 )
{
force(a, y, nB);
h = dt * d8[i];
for(int j = 0; j < nB; j++)
y[j].v.addm(a[j], h);
}
}
}

but the extrapolated method - even limited to the 8-order, is still much better than this.