C++でリープフロッグ法

リープフロッグ(LeapFrog)法のプログラム。
ケプラー問題

を解いた。
G=1、M=1として初期条件をr(1, 0, 0)、v(0, 1, 0)とした。

#include<iostream>
#include<math.h>
 
using namespace std;
 
#define VLEN 3
#define REAL double
 
class vector{
private:
  REAL element[VLEN];
public:
  vector(REAL c = 0){
    for(int i = 0; i<VLEN; i++){
      element[i] = c;
    }
  }
  double & operator [] (int i){
    return element[i];
  }
 
  vector & operator += (const vector & b){
    for(int i = 0; i < VLEN; i++){
      element [i] += b.element[i];
    }
    return *this;
  }
 
  friend const vector operator + (const vector &, const vector &);
  friend const vector operator - (const vector &, const vector &);
  friend const vector operator * (const REAL &, const vector &);
  friend const double operator * (const vector &, const vector &);
  friend const vector operator | (const vector &, const vector &);
  friend ostream & operator << (ostream &, const vector &);
  friend istream & operator >> (istream &, vector &);
 
};
 
inline ostream & operator << (ostream & s, const vector & v){
  for(int i = 0; i < VLEN; i++){
    s << v.element[i] << " ";
  }
  return s;
}
 
inline istream & operator >> (istream & s, vector & v){
  for (int i = 0; i < VLEN; i++){
    s >> v.element[i];
  }
  return s;
}
 
inline const vector operator + (const vector &v1, const vector &v2){
  vector v3;
  for(int i = 0; i < VLEN; i++ ){
    v3.element[i] = v1.element[i] + v2.element[i];
  }
  return v3;
}
 
inline const vector operator - (const vector &v1, const vector &v2){
  vector v3;
  for(int i = 0; i < VLEN; i++ ){
    v3.element[i] = v1.element[i] - v2.element[i];
  }
  return v3;
}
 
inline const vector operator * (const REAL &c, const vector &v1){
  vector v3;
  for(int i = 0; i < VLEN; i++ ){
    v3.element[i] = v1.element[i] * c;
  }
  return v3;
}
 
inline const double operator * (const vector &v1, const vector &v2){
  double sum = 0;
  for(int i = 0; i < VLEN; i++ ){
    sum += v1.element[i] * v2.element[i];
  }
  return sum;
}
 
inline const vector operator | (const vector &v1, const vector &v2){
  vector v3;
  v3[0] = v1.element[1] * v2.element[2] - v1.element[2] * v2.element[1];
  v3[1] = v1.element[2] * v2.element[0] - v1.element[0] * v2.element[2];
  v3[2] = v1.element[0] * v2.element[1] - v1.element[1] * v2.element[0];
  return v3;
}
 
void leapfrog(vector&, vector&, vector&, REAL);
void func(vector, vector&);
 
main(){
  FILE * data;
  int n;
  vector r, v, e, enew, de, u;
  REAL t, tend, h, nstep, l, a, energy, lnew, anew, energynew;
  REAL da, denergy, dl, p;
 
  r[0] = 1;
  r[1] = 0;
  r[2] = 0;
  v[0] = 0;
  v[1] = 1;
  v[2] = 0;
 
  h = 0.1;
  t = 0.0;
  data = fopen("4.2-leapfrog.dat", "w");
 
  e = (v|(r|v)) - ((1.0/(sqrt(r*r)))*r);
  energy = (1.0/2.0)*v*v -(1.0/(sqrt(r*r)));
  a = -1.0/(2.0*energy);
  l = sqrt(a*(1.0-e*e));
 
  p = 2*M_PI*a*a/sqrt(a);
  tend = 10 * p;
  nstep =1 + (int)(tend/h);
 
  fprintf(data, "%8.3f %16.8e %16.8e %16.8e %16.8e ", t,energy,l,e*e,a);
  fprintf(data, "%16.8e %16.8e %16.8e " ,r[0],r[1],r[2]);
  fprintf(data, "%16.8e %16.8e %16.8e \n", e[0],e[1],e[2]);
 
  func(r, u);  
  for(int i = 0; i < nstep; i++){
    leapfrog(r, v, u, h);
    t += h;
    if(i % 10 == 9){
      enew = (v|(r|v)) - ((1/(sqrt(r*r)))*r);
      energynew = (1.0/2.0)*v*v -(1.0/(sqrt(r*r)));
      anew = -1.0/(2.0*denergy);
      lnew = sqrt(da*(1.0-enew*enew));
 
      de = enew - e;
      denergy = energynew -energy;
      da = anew - a;
      dl = lnew - l;
 
      fprintf(data, "%8.3f %16.8e %16.8e %16.8e %16.8e ", t,denergy,dl,de*de,da);
      fprintf(data, "%16.8e %16.8e %16.8e ", r[0],r[1],r[2]);
      fprintf(data, "%16.8e %16.8e %16.8e \n", de[0],de[1],de[2]);
    }
  }
  fclose(data);
}
 
void leapfrog(vector &r, vector &v, vector &u, REAL h){
  static vector v2;
  v2 = v + 0.5*h*u;
  r = r + h*v2;
  func(r, u);
  v = v2 + 0.5*h*u;
}
 
void func(vector x, vector &u1){
  double a = 1;
  double b;
  b = x * x;
  u1 = ((-1/sqrt(b))/b) * x;
}
 

20周期でやった結果。