Dismiss Notice
Join Physics Forums Today!
The friendliest, high quality science and math community on the planet! Everyone who loves science is here!

C/++/# Runge-Kutta simulation in C++

  1. Mar 31, 2016 #1
    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)

    Code (C):

    #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>
     

    Attached Files:

    Last edited by a moderator: Mar 31, 2016
  2. jcsd
  3. Mar 31, 2016 #2

    ShayanJ

    User Avatar
    Gold Member

    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.
     
  4. Mar 31, 2016 #3
    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?
     
  5. Mar 31, 2016 #4

    DrClaude

    User Avatar

    Staff: Mentor

    Do you know what the 4 stands for in RK4?
     
  6. Mar 31, 2016 #5

    ShayanJ

    User Avatar
    Gold Member

    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!
     
  7. Mar 31, 2016 #6

    Filip Larsen

    User Avatar
    Gold Member

    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
     
  8. Mar 31, 2016 #7
    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
     
  9. Mar 31, 2016 #8

    jtbell

    User Avatar

    Staff: Mentor

    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.
     
  10. Mar 31, 2016 #9
    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?
     
  11. Apr 2, 2016 #10

    D H

    User Avatar
    Staff Emeritus
    Science Advisor

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

    Accuracy and precision are distinct issues.
    accuracy_vs_precision_556.jpg

    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?
     
  12. Apr 3, 2016 #11
    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).
     
  13. Apr 4, 2016 #12
    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.
     
  14. Apr 4, 2016 #13

    Mark44

    Staff: Mentor

    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.
     
  15. Apr 4, 2016 #14

    jim mcnamara

    User Avatar

    Staff: Mentor

    @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.
     
  16. Apr 5, 2016 #15
    There are also lots of libraries that define __float128
     
Know someone interested in this topic? Share this thread via Reddit, Google+, Twitter, or Facebook

Have something to add?
Draft saved Draft deleted