N-body仿真不能正常工作

Question

我正在尝试编写一个代码来解决N-body问题，并且在与所有机构一起使用数组而不是单独使用不同的机构时遇到了麻烦。我目前不知道我的代码中出了什么问题，但是当我在任何物体的y函数中绘制x时，我会得到一条明显不正确的直线。

这是我当前的代码：

#include <cstdlib> 
#include <iostream> 
#include <cmath> 
#include <fstream> 
#define h 10000.0 
#define N 3 
#define G 6.67384*pow(10.0,-11) 

using namespace std; 



class particle{ 
     public: 
     double kx1,kx2,kx3,kx4, kv1, kv2, kv3, kv4; 
     double ky1, ky2, ky3, ky4, kvy1, kvy2, kvy3, kvy4; 
     double x,y,vx,vy,m; 


     double dist(particle aap){ 
      double dx = x - aap.x; 
      double dy = y - aap.y; 
      return sqrt(pow(dx,2.0)+pow(dy,2.0)); 
      } 

     double g(double x1, double y1,particle aap){ 
      return G*aap.m*(aap.x-x1)/pow(dist(aap),3.0); 
      } 

     double p(double x1, double y1, particle aap){ 
      return G*aap.m*(aap.y-y1)/pow(dist(aap),3.0); 
     } 


     void update(){          //object advances 1 step 
      x = x + (1/6.0)*(kx1+2*kx2+2*kx3+kx4); 
      vx = vx + (1/6.0)*(kv1+2*kv2+2*kv3+kv4); 
      y = y + (1/6.0)*(ky1+2*ky2+2*ky3+ky4); 
      vy = vy + (1/6.0)*(kvy1+2*kvy2+2*kvy3+kvy4); 
      } 

    void create(double x1, double y1, double vx1, double vy1, double m1){ 
         x = x1; 
         y = y1; 
         vx = vx1; 
         vy = vy1; 
         m =m1; 
         } 

    bool operator ==(particle &other){ 
      if(x == other.x && y == other.y && vx == other.vx && vy == other.vy){ 
       return true; 
       } 
       } 

     }; 




particle bodies[N]; 



void set(particle (&bodies)[N]){ 
    bodies[0].create(1, 1, -2, 1, 2*pow(10.0,30)); 
    bodies[1].create(2870671*pow(10.0,6), 0, 0, 6800, 8.6810*pow(10.0,25)); 
    bodies[2].create(4498542*pow(10.0,6),0 ,0, 5430, 1.0243*pow(10.0,26)); 
    } 


double xforce(double x1, double y1, particle aap, particle bodies[N]){ //force in the x- direction 

     double fx = 0; 
     for (int i = 0; i <= N; i++){ 
      if (bodies[i] == aap){;} 

      else{ 
       fx += aap.g(x1,y1,bodies[i]); 
       } 
       } 
     return fx; 
     } 

double yforce(double x1, double y1, particle aap, particle bodies[N]){ //force in the y- direction 

     double fy = 0; 
     for (int i = 0; i <= N; i++){ 
      if (bodies[i] == aap) {;} 

      else{ 
       fy += aap.p(x1,y1,bodies[i]); 
       } 
       } 
     return fy; 
     } 


void corr(double t, particle bodies[N]){       //runge kutta 4 
    for(int i =0; i <= N; i++){ 

      bodies[i].kx1 = t*bodies[i].vx; 
      bodies[i].kv1 = t*xforce(bodies[i].x, bodies[i].y, bodies[i], bodies); 
      bodies[i].ky1 = t*bodies[i].vy; 
      bodies[i].kvy1 = t*yforce(bodies[i].x, bodies[i].y, bodies[i], bodies); 

      bodies[i].kx2 = t*(bodies[i].vx + 0.5*bodies[i].kv1); 
      bodies[i].kv2 = t*xforce(bodies[i].x + 0.5*bodies[i].kx1, bodies[i].y + 0.5*bodies[i].ky1, bodies[i], bodies); 
      bodies[i].ky2 = t*(bodies[i].vy + 0.5*bodies[i].kvy1); 
      bodies[i].kvy2 = t*yforce(bodies[i].x + 0.5*bodies[i].kx1, bodies[i].y + 0.5*bodies[i].ky1, bodies[i], bodies); 

      bodies[i].kx3 = t*(bodies[i].vx+ 0.5*bodies[i].kv2); 
      bodies[i].kv3 = t*xforce(bodies[i].x + 0.5*bodies[i].kx2, bodies[i].y + 0.5*bodies[i].ky2, bodies[i], bodies); 
      bodies[i].ky3 = t*(bodies[i].vy+ 0.5*bodies[i].kvy2); 
      bodies[i].kvy3 = t*yforce(bodies[i].x + 0.5*bodies[i].kx2, bodies[i].y + 0.5*bodies[i].ky2,bodies[i], bodies); 

      bodies[i].kx4 = t*(bodies[i].vx + bodies[i].kv3); 
      bodies[i].kv4 = t*xforce(bodies[i].x+ bodies[i].kx3, bodies[i].y + bodies[i].ky3, bodies[i], bodies); 
      bodies[i].ky4 = t*(bodies[i].vy + bodies[i].kvy3); 
      bodies[i].kvy4 = t*yforce(bodies[i].x + bodies[i].kx3, bodies[i].y + bodies[i].ky3, bodies[i], bodies); 
      } 
    } 


void calculate(particle (&bodies)[N]){ 
    set(bodies); 
    ofstream file; 
    file.open("tester.txt"); 
    for(int i =0; i <=50000; i++){ 

      corr(h, bodies);       
      for(int j = 0; j <= N; j++){ 
        bodies[j].update(); 
        }     
      if(i%1000 == 0){ 

       file << i*h; 
       for(int j = 0; j <=N ; j++){ 
          file <<" "<<bodies[j].x << " "<< bodies[j].y; 
          } 
       file <<" "<<"\n"; 
       } 
      else{;} 
      } 
    file.close(); 
    } 


int main() 
{ 
    calculate(bodies); 
    system("pause"); 
    return 0; 
} 

的问题可能存在的类颗粒之外，因为该计划的工作之前，我开始使用数组机构。任何关于非必要改进的建议都是可喜的。我想要做的另一件事是使用std :: vector而不是数组，但我不知道如何定义一个向量之外的向量，就像我定义的数组体。

Answer 1

首先，您所有的i <= N都是错误的，因为您的循环将执行4次（0,1,2,3）而不是3次，用于i < N。