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.
Pages: 1 2
faster and better method of integration (Read 9798 times)
