#include
#include
#include
/*
* symp.cpp
* Symplectic integration routines to 4th, 6th, and 8th order.
* Copyright (c) 1998 David Whysong, with alterations by Bill Gray.
* (Further revised 2013 April 26 by Bill Gray: comments added,
* 'const' declarations added, more digits shown in results, revised
* to compile cleanly in Linux.)
*
* This program is free software; you can redistribute it and/or modify
* it under the terms of the GNU General Public License as published by
* the Free Software Foundation, version 2.
*
* Reference: "Construction of higher order symplectic integrators" Haruo
* Yoshida, Phys. Lett. A 150, pp. 262, 12 Nov. 1990.
* http://cacs.usc.edu/education/phys516/Yoshida-symplectic-PLA00.pdf
*
* For the test program, compile as:
*
* cl -W4 -Ox -DTEST_CODE symp.cpp (Microsoft C/C++)
* wcl386 -W4 -Ox -DTEST_CODE symp.cpp (Watcom C/C++)
* g++ -Wall -o symp -O3 -DTEST_CODE symp.cpp
*
* The test case comes from p. 168, J M A Danby's _Fundamentals of
* Celestial Mechanics_; it corresponds to an object making about
* one-sixth of a roughly circular orbit. The test code shows the
* ending state vector, plus the error (numerical solution minus
* analytical solution). You run the test program with a number of
* steps and a '4', '6', or '8' to indicate the method to be used;
* for example,
*
* symp 19 6
*
* to use the symplectic_6 routine, with 19 integration steps. You'll
* notice that doubling the number of steps cuts the error for
* symplectic_4 about 16-fold; for symplectic_6, 64-fold; and for
* symplectic_8, about 256-fold (which matches the fact that these
* routines should be fourth, sixth, and eighth-order, respectively.)
*
* NOTE that at first blush, the idea of using a sixth or eighth-order
* integrator is going to look wonderful. However, the fourth-order method
* uses three force-function evaluations; the sixth-order uses seven; and
* the eighth-order uses 15. Thus, for example, all of these :
*
* ./symp 28 8
* ./symp 140 4
* ./symp 60 6
*
* should require 28*15 = 140*3 = 60*7 = 420 force-function evaluations,
* and will therefore take about the same amount of computational time.
* At least in this case, the sixth-order method provides smaller errors
* than the other two. But if you needed higher accuracy (and therefore
* increased the number of steps), eventually, the eighth-order method
* would be the winner; and for a sufficiently rough evaluation, you
* would find that the fourth-order method worked best.
*
* The 'acceleration' code takes the time and velocity as parameters,
* but in this particular case, the acceleration is a function of
* position only. I'm pretty sure I've got the time and velocity
* dependence right, but can't say for an absolute certainty. If
* one wanted to include planetary perturbations (which vary with
* time) or some sort of drag term (a function of velocity), this
* would suddenly matter.
*
* I'll throw in some comments on why this is useful stuff, derived
* mostly from e-mail talks with David Whysong. The big plus with a
* symplectic method is that quantities such as angular momentum and
* total energy are pretty well-conserved (more so than with other
* methods such as Runge-Kutta, Adams-Bashforth, Adams-Bashforth-
* Moulton, and their kin.) David says they have only really emerged
* in the last decade or so, which explains why I never heard of them
* in physics courses in 1983-87. He's using them in research of
* globular cluster simulations, where it's necessary to run _very_
* long integrations. If, say, energy was not conserved over so long
* an integration, accumulated errors might cause your simulated
* globular cluster to undergo a simulated explosion, or a simulated
* gravitational collapse. I _don't_ know if a price is paid in
* precision for this benefit. I also have found no references to
* symplectic integration.
*
*/
void compute_accel( double *accel, const double *x, const double *v,
const double t)
{
int i;
const double r2 = x[0] * x[0] + x[1] * x[1] + x[2] * x[2];
const double a0 = 1. / (r2 * sqrt( r2));
for( i = 0; i < 3; i++)
accel[i] = -x[i] * a0;
}
void symplectic_4( double *x, double *v, double t, const double dt)
{
int i, j;
const double b = 1.25992104989487319066654436028; /* 2^1/3 */
const double a = 2 - b;
const double x0 = -b / a;
const double x1 = 1. / a;
const static double d4[3] = { x1, x0, x1 };
const static double c4[4] = { x1/2, (x0+x1)/2, (x0+x1)/2, x1/2 };
for( i = 0; i < 4; i++)
{
double accel[3];
const double step = dt * c4[i];
for( j = 0; j < 3; j++)
x[j] += step * v[j];
t += step;
if( i != 3)
{
compute_accel( accel, x, v, t);
for( j = 0; j < 3; j++)
v[j] += dt * d4[i] * accel[j];
}
}
}
void symplectic_6( double *x, double *v, double t, const double dt)
{
int i, j;
const double w1 = -0.117767998417887E1;
const double w2 = 0.235573213359357E0;
const double w3 = 0.784513610477560E0;
const double w0 = (1-2*(w1+w2+w3));
const static double d6[7] = { w3, w2, w1, w0, w1, w2, w3 };
const static double c6[8] = { w3/2, (w3+w2)/2, (w2+w1)/2, (w1+w0)/2,
(w1+w0)/2, (w2+w1)/2, (w3+w2)/2, w3/2 };
for( i = 0; i < 8; i++)
{
double accel[3];
const double step = dt * c6[i];
for( j = 0; j < 3; j++)
x[j] += step * v[j];
t += step;
if( i != 7)
{
compute_accel( accel, x, v, t);
for( j = 0; j < 3; j++)
v[j] += dt * d6[i] * accel[j];
}
}
}
void symplectic_8( double *x, double *v, double t, const double dt)
{
int i, j;
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( i = 0; i < 16; i++)
{
double accel[3];
const double step = dt * c8[i];
for( j = 0; j < 3; j++)
x[j] += step * v[j];
t += step;
if( i != 15)
{
compute_accel( accel, x, v, t);
for( j = 0; j < 3; j++)
v[j] += dt * d8[i] * accel[j];
}
}
}
#ifdef TEST_CODE
int main( const int argc, const char **argv)
{
double x[6] = {1., .1, -.1, -.1, 1., .2};
const double analytic[6] = {
.4660807846539234, .9006112349077236, .1140459806673234,
-.8765944368525084, .4731566049794786, .1931594924439623 };
const int n_steps = atoi( argv[1]);
int i;
for( i = 0; i < n_steps; i++)
switch( argv[2][0])
{
case '4':
symplectic_4( x, x + 3, 0., 1. / (double)n_steps);
break;
case '6':
symplectic_6( x, x + 3, 0., 1. / (double)n_steps);
break;
case '8':
symplectic_8( x, x + 3, 0., 1. / (double)n_steps);
break;
default:
printf( "Option '%s' not understood\n", argv[2]);
return( -1);
break;
}
printf( "position: %18.15lf %18.15lf %18.15lf \n", x[0], x[1], x[2]);
printf( "velocity: %18.15lf %18.15lf %18.15lf \n", x[3], x[4], x[5]);
printf( "posn err: %18.15lf %18.15lf %18.15lf \n",
x[0] - analytic[0], x[1] - analytic[1], x[2] - analytic[2]);
printf( "vel err: %18.15lf %18.15lf %18.15lf \n",
x[3] - analytic[3], x[4] - analytic[4], x[5] - analytic[5]);
return( 0);
}
#endif