www.pudn.com > symplectic.rar > symplectic.cpp, change:2006-08-05,size:6083b


/* 
 *  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 <dwhysong@physics.ucsb.edu> 
 * 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