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。