Why Do My Simulated Orbits Appear as Spirals in C++ Runge-Kutta Simulation?

[Silviu]

Hi! I have this code in C++ that simulates 2 body interacting through gravity using Runge-Kutta 4th order method. I plotted the trajectory of the 2 bodies and as you can see in the first picture they seems to be what I expect, 2 circles. But when I zoom in on one of them (second image) it looks like the circle is not really a circle but more like a spiral. Could you please tell me what might be wrong? Thank you! (The plot is done with python)

#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,ax,ay;

    double fx;

    double fy;

    double fz;

    double dt = 0.1;

    double initial_x, initial_y,initial_vx,initial_vy;

    double dist,dx,dy,dz;

    double x1,x2,x3,x4,vx1,vx2,vx3,vx4,ax1,ax2,ax3,ax4;

    double y1,y2,y3,y4,vy1,vy2,vy3,vy4,ay1,ay2,ay3,ay4;

  

  

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){

        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);

      

        vx2=vx+ax1*dt*0.5;

        vy2=vy+ay1*dt*0.5;

      

        x2=rx+vx*dt*0.5;

        y2=ry+vy*dt*0.5;

      

        dx=b.rx-x2;

        dy=b.ry-y2;

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

      

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

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

      

        vx3=vx+ax2*dt*0.5;

        vy3=vy+ay2*dt*0.5;

      

        x3=rx+vx2*dt*0.5;

        y3=ry+vy2*dt*0.5;

      

        dx=b.rx-x3;

        dy=b.ry-y3;

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

      

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

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

      

        vx4=vx+ax3*dt;

        vy4=vy+ay3*dt;

      

        x4=rx+vx3*dt;

        y4=ry+vy3*dt;

      

        dx=b.rx-x4;

        dy=b.ry-y4;

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

      

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

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

      

        double nr=1.0/6;

      

        rx=rx+(nr)*(vx+2*vx2+2*vx3+vx4)*dt;

        ry=ry+(nr)*(vy+2*vy2+2*vy3+vy4)*dt;

      

        vx=vx+(nr)*(ax1+2*ax2+2*ax3+ax4)*dt;

        vy=vy+(nr)*(ay1+2*ay2+2*ay3+ay4)*dt;

      

      

    }

  

  

  

    double get_x(){

        return rx;

    }

  

    double get_y(){

        return ry;

    }

  

    double get_vx(){

        return vx;

    }

  

};

int main(){
    Body body1(0,0,0,0,-0.527046,0,1000), body2(600,0,0,0,1.05409,0,500);

  

    ofstream pos;

    pos.open ("Position.txt");

  

  

    int N=100000;

    for(int i; i<N;i++){

        body2.update(body1);
        body1.update(body2);

        pos<<body2.get_x()<<" "<<body2.get_y()<<" "<<body1.get_x()<<" "<<body1.get_y()<<endl;

    }

  

  

    pos.close();

  

}

<Moderator's note: Please use code tags when posting code>

 

[ShayanJ]

You always have a finite precision which in this case is controlled by your choice of method and the value of dt. RK4 is good enough for this problem so just keep playing with dt. The smaller it gets, the more precise will become your plot.

 

[Silviu]

You always have a finite precision which in this case is controlled by your choice of method and the value of dt. RK4 is good enough for this problem so just keep playing with dt. The smaller it gets, the more precise will become your plot.

Thank you! But I need to simulate numerous periods (the system will be more complex than 2 bodies) so even if I make dt smaller I will have to increase the number of iterations to get the proper number of rotations. So in the end won't the result be wrong even with a very small dt?

 

[DrClaude]

Thank you! But I need to simulate numerous periods (the system will be more complex than 2 bodies) so even if I make dt smaller I will have to increase the number of iterations to get the proper number of rotations. So in the end won't the result be wrong even with a very small dt?

Do you know what the 4 stands for in RK4?

 

[ShayanJ]

Thank you! But I need to simulate numerous periods (the system will be more complex than 2 bodies) so even if I make dt smaller I will have to increase the number of iterations to get the proper number of rotations. So in the end won't the result be wrong even with a very small dt?

The result is not wrong! You're looking at it the wrong way.
As I said before, you always have a finite precision in such numerical computations. You should decide to what precision you want to get the answer, and write your code somehow that you get that precision.
And about the accumulation error, I guess its better I leave it to DrClaude to explain it!

 

[Filip Larsen]

The effect you observe is consistent with energy drift [1] which will be present in the numerical solution for Hamiltonian systems when the integrator used is not symplectic [2]. Note, that no explicit Runge-Kutta method is symplectic so if you want to avoid energy drift you may want to investigate methods like Verlet integration [3] or other symplectic integrators.

[1] https://en.wikipedia.org/wiki/Energy_drift
[2] https://en.wikipedia.org/wiki/Symplectic_integrator
[3] https://en.wikipedia.org/wiki/Verlet_integration

 

[Silviu]

Do you know what the 4 stands for in RK4?

I am new to this, but from what I understood is related to the error propagation (like it goes like dt^4 or something). Is it right? And thank you for trying to help

 

[jtbell]

I understood is related to the error propagation (like it goes like dt^4 or something)

Basically, yes. As you may be aware, there is a 2nd order Runge-Kutta method which is simpler but less accurate. There are also higher-order Runge-Kutta methods which are more accurate. However, they will simply reduce the systematic error that you see, not eliminate it. To eliminate this particular kind of systematic error, you need to switch to a different kind of method e.g. Verlet, as others have mentioned.

 

[Silviu]

Basically, yes. As you may be aware, there is a 2nd order Runge-Kutta method which is simpler but less accurate. There are also higher-order Runge-Kutta methods which are more accurate. However, they will simply reduce the systematic error that you see, not eliminate it. To eliminate this particular kind of systematic error, you need to switch to a different kind of method e.g. Verlet, as others have mentioned.

Thank you! Two more question, how accurate is Verlet compared to Runge-Kutta? (I want a precise result). And what are the advantages of using Runge-Kutta compared to Verlet, if Verlet doesn't give that kind of errors?

 

[D H]

And what are the advantages of using Runge-Kutta compared to Verlet, if Verlet doesn't give that kind of errors?

Accuracy, precision, and significantly less computation, at least over short spans of time, where short means several dozen orbits.

Two more question, how accurate is Verlet compared to Runge-Kutta? (I want a precise result).

Accuracy and precision are distinct issues.

With regard to numerically integrating a second order differential equation such as gravitation, omputational effort is yet another issue, code complexity, yet another, and stability, yet another. You can't maximize four of them, let alone all five. Addressing three is not easy. So what do you want?

 

[JorisL]

Even with Verlet integration you'll see various lines if you zoom in.
In fact the computed position oscillate around the true orbit, in principle this is true no matter how long you simulate (if I'm not mistaken).

 

[newjerseyrunner]

Replace your "double" with "long double"
Don't forget to change your constants:
double G=1;
should be
long double G = 1L;

It will add precision to your calculations, but as others have said, it will not be perfect.

 

[Mark44]

Replace your "double" with "long double"

That may or may not help. If the OP is using the Microsoft C/C++ compiler, long double and double are the same size, eight bytes. The GNU C compiler has a long double type of 80 bits that is larger than the double type.

 

[jim mcnamara]

@Silviu

I think you have a wrong assumption about arithmetic operations with floating point numbers. Instead of thrashing about consider this as one of the things you need to know:

https://www.physicsforums.com/insights/cant-computer-simple-arithmetic/

Floating point simply cannot accurate represent some numbers - for example, if you were working on balancing ledgers in bank (money) you would NOT want to use double or long double. Auditors will call you on it. And then you get to change all of your code...

Also IMO, going to long double simply moves the fuzz further out. If you zoom in even further you will find it is still there.

PS: what you did is NOT wrong, it just has very, very slightly fuzzy numbers. In some cases.

 

[newjerseyrunner]

That may or may not help. If the OP is using the Microsoft C/C++ compiler, long double and double are the same size, eight bytes. The GNU C compiler has a long double type of 80 bits that is larger than the double type.

There are also lots of libraries that define __float128