symplectic.cpp

4阶6阶及8阶的辛几何算法源程序计算分子模型

/*
 *  symplectic.cpp
 *  //Symplectic integration routines to 4th, 6th, and 8th order.//
 *  Copyright (C) 1998, 1999 David Whysong with alterations by Bill Gray.
 *
 *   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.
 *
 * David Whysong <[e-mail]>
 * March 4 1999
 * http://www.physics.ucsb.edu/~dwhysong
 *
 *  For the test program,  compile as:
 *
 *    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'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.)
 *
 *    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.
 *
 */

#include <stdio.h>
#include <stdlib.h>
#include <math.h>

void compute_accel(double *accel, double *x, double *v, double t)
{
   int i;
   double r2 = x[0] * x[0] + x[1] * x[1] + x[2] * x[2];
   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, 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];
      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, 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];
      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, 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];
      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];
         }
      }
}

#ifndef TEST_CODE
#define TEST_CODE

int main(int argc, char **argv)
{
   double x[6] = {1., .1, -.1, -.1, 1., .2};
   double analytic[6] = {
          .4660807846539234, .9006112349077236, .1140459806673234,
         -.8765944368525084, .4731566049794786, .1931594924439623 };
   int n_steps;
   int i;
   int order;
	printf("Enter an integer steps >" );
	scanf("%d",& n_steps);
	printf("Enter an integer order of integration  >\n" 
			 "Available integrators are 4th, 6th, and 8th order.\n\n");
	scanf("%d",&order);


   for( i = 0; i < n_steps; i++)
     	  switch( order)
         {
         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;
         }
   printf("position:  %15.12g %15.12g %15.12g \n", x[0], x[1], x[2]);
   printf("velocity:  %15.12g %15.12g %15.12g \n", x[3], x[4], x[5]);
   printf("posn err:  %15.12g %15.12g %15.12g \n",
         x[0] - analytic[0], x[1] - analytic[1], x[2] - analytic[2]);
   printf("vel err:   %15.12g %15.12g %15.12g \n",
         x[3] - analytic[3], x[4] - analytic[4], x[5] - analytic[5]);
   return(0);
}
#endif