C/C++ Runge-Kutta methond for n body C++

[Silviu]

Hi! I am trying to simulate the rotation of a planet around a star, using the 2nd order Runge-Kutta method (I am starting with this and I will try 4th order later, I am new to this topic). I have this code but it doesn't work. When I plot y(x) I should get a circle, but I don't get it. Any idea what am I doing wrong? Thank you!

#include<iostream>
#include <vector>
#include<math.h>
#include <fstream>

usingnamespace std;

class Body{
  private:
    double G= 1;
    double rx;
    double ry;
    double rz;
    double vx;
    double vy;
    double vz;
    double mass;
    double dt = 0.001;
    double dist,dx,dy;
    double x1,vx1,ax1,ax2;
    double y1,vy1,ay1,ay2;

  public:
    Body(double rx, double ry, double rz, double vx, double vy, double vz, double mass){
        this->rx=rx;
        this->ry=ry;
        this->rz=rz;
        this->vx=vx;
        this->vy=vy;
        this->vz=vz;
        this->mass=mass;
    }

    void update(Body b){

        //Step 1
        dx=b.rx-rx;
        dy=b.ry-ry;
        dist = sqrt(dx*dx+dy*dy);

        ax1=G*b.mass*dx/(dist*dist*dist);
        ay1=G*b.mass*dy/(dist*dist*dist);
  
        //Step 2
        x1=rx+vx*dt*0.5;
        vx1=vx+ax1*dt*0.5;

        y1=ry+vy*dt*0.5;
        vy1=vy+ay1*dt*0.5;

        dx=b.rx-x1;
        dy=b.ry-y1;

        dist = sqrt(dx*dx+dy*dy);

        ax2=G*b.mass*dx/(dist*dist*dist);
        ay2=G*b.mass*dy/(dist*dist*dist);

        rx=rx+0.5*(vx1+vx)*dt;
        ry=ry+0.5*(vy1+vy)*dt;

        vx=vx+0.5*(ax1+ax2)*dt;
        vy=vy+0.5*(ay1+ay2)*dt;   
    }

    double get_x(){
        return rx;
    }

    double get_y(){
        return ry;
    }
};

int main(){
    Body body1(0,0,0,0,0,0,1000), body2(10,0,0,0,10,0,10);
    ofstream pos;

    pos.open ("Position.txt");
    int N=100000;

    for(int i; i<N;i++){
        body2.update(body1);
        pos<<body2.get_x()<<" "<<body2.get_y()<<endl;
    }

    pos.close();
}

<<Mentor's note: please use code tags when posting code.>>

 

[Twigg]

I spotted at least one major typo that would be throwing your simulation off:

ay2=G*b.mass*dx/(dist*dist*dist);

Change that to:

ay2=G*b.mass*dy/(dist*dist*dist);

You should look for additional typos.

More importantly, since you say you're new to the topic, you should take notes how troubleshooting this code took a long time to resolve a simple issue. One of the essential skills of a programmer in any field is writing code that's easy to fix. The most important technique is proper commenting. You don't need a billion comments, you just need enough. I would advise you to get into the habit of writing comments before every method declaration that contain at least one example of how the method works. For example:

// The redundant(a,b) function takes a character input a and prints it to the screen some number of times b.
// Inputs: a--some character to print
// b--an integer value, how many times to print the character a
// Example: redundant('a',3) prints "aaa" to standard output

void redundant(char a,int b) {
for (int ii=0;ii<b;ii++) {
std::cout<<a;
}
}

It looks tedious, and sometimes it is, but it's still less tedious than debugging uncommented code. There are other organizational tricks I would recommend to make your life easier, but if you can just get this one habit down you'll be a happier coder. Especially if you want to do something as intricate as many body problems. I hope you find this helpful.

 

[Lord Crc]

I would also suggest you use a simple problem which has a simple analytical solution first, if you haven't already.

For example the 1D ballistic problem is nice because you can very easily verify that the velocity and acceleration makes sense for each step, by printing the actual values and the expected values.