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

Two Body problem with RK4 in C++

  1. Dec 14, 2012 #1
    I want to output an excel file with the results of the trajectories of a two body problem, with initial position and velocity. But my program is not compiling. Any suggestions/problems that you can see?

    Code (Text):

    #include <iostream>
    #include <vector>
    #include <fstream>

    using namespace std;

    void RK4_func(double,double,double,double,double,double,double,double,double,double,double,double);

    int main(){


    //initialise variables
        double time_end,time_step,mass1,x1x,y1y,v1x,v1y,mass2,x2x,y2y,v2x,v2y;

    //Get the value
        cout << "Enter Simulation Duration in seconds" << endl;
        cin >> time_step;
        cout << "\n";
       
        cout << "Enter Precision (Time step)" << endl;
        cin >> time_end;
        cout << "\n";
       
        cout << "Enter mass of body 1 with 10^20kg scale" << endl;
        cin >> mass1;
        cout << "\n";
       
        cout << "Enter initial body 1 x position" << endl;
        cin >> x1x;
        cout << "\n";
       
        cout << "Enter initial body 1 y position" << endl;
        cin >> y1y;
        cout << "\n";
       
        cout << "Enter initial body 1 x velocity" << endl;
        cin >> v1x;
        cout << "\n";
       
        cout << "Enter initial body 1 x velocity" << endl;
        cin >> v1y;
        cout << "\n";
       
        cout << "Enter mass of body 2 with 10^20kg scale" << endl;
        cin >> mass2;
        cout << "\n";
       
        cout << "Enter initial body 2 x position" << endl;
        cin >> x2x;
        cout << "\n";
       
        cout << "Enter initial body 2 y position" << endl;
        cin >> y2y;
        cout << "\n";
       
        cout << "Enter initial body 2 x velocity" << endl;
        cin >> v2x;
        cout << "\n";
       
        cout << "Enter initial body 2 y velocity" << endl;  
        cin >> v2y;
        cout << "\n";
       

    //Pass values int RK4_func
    mass1=mass1*10^20
    mass2=mass2*10^20
    RK4_func(time_end,time_step;mass1;x1x,y1y,v1x,v1y,mass2,x2x,y2y,v2x,v2y);

    //Output data

        return 0;
    }

    //Assume function  enters data in SI units

    void RK4_func (double timeend, double steps, double m1, double x1, double y1, double vx1, double vy1, double m2, double x2, double y2, double vx2,double vy2)
    {

    //Excel file creation

        ofstream outfile1 ("RK4.xls"); 
    outfile1<<"time"<<'\t'<<"x1"<<'\t'<<"y1"<<'\t'<<"vx1"<<'\t'<<"vy1"<<'\t'<<"x2"<<'\t'<<"y2"<<'\t'<<"vx2"<<'\t'<<"vy2"<<endl;
    outfile1<<0<<'\t'<<x1<<'\t'<<y1<<'\t'<<vx1<<'\t'<<vy1<<'\t'<<x2<<'\t'<<y2<<'\t'<<vx2<<'\t'<<vy2<<endl;

        constant double G = 6.67E-11;
        double datalength = 1+(timeend/steps);
        double data [8][datalength];
       
        double t=0;
        data[0][0]=x1;
        data[1][0]=y1;
        data[2][0]=x2;
        data[3][0]=y2;
        data[4][0]=vx1;
        data[5][0]=vy1;
        data[6][0]=vx2;
        data[7][0]=vy2;

        for (int i=0; i<datalength; i++)
        {
            vector<double> interdata(12), k1 (8), k2 (8), k3 (8), k4 (8);
       
            // let 0~3 indicate position x1y1x2y2 4~7 indicate vx1vy1vx2vy2 8~11 indicate ax1ay1ax2ay2
           
            //acceleration x1 y1 x2 y2
           
            interdata[0]=data[0][i];
            interdata[1]=data[1][i];
            interdata[2]=data[2][i];
            interdata[3]=data[3][i];
            interdata[4]=data[4][i];
            interdata[5]=data[5][i];
            interdata[6]=data[6][i];
            interdata[7]=data[7][i];
            interdata[8] = G*m2*/(((((data[0][i]-data[2][i])^2)+((data[1][i]-data[3][i])^2))^(1/2))^3)*(data[2][i]-data[0][i]);
            interdata[9] = G*m2*/(((((data[0][i]-data[2][i])^2)+((data[1][i]-data[3][i])^2))^(1/2))^3)*(data[3][i]-data[1][i]);
            interdata[10] = G*m1*/(((((data[0][i]-data[2][i])^2)+((data[1][i]-data[3][i])^2))^(1/2))^3)*(data[0][i]-data[2][i]);
            interdata[11] = G*m1*/(((((data[0][i]-data[2][i])^2)+((data[1][i]-data[3][i])^2))^(1/2))^3)*(data[1][i]-data[3][i]);
           
            //k1 is Euler Step
           
            k1[0] = steps*interdata[8];
            k1[1] = steps*interdata[9];
            k1[2] = steps*interdata[10];
            k1[3] = steps*interdata[11];
            k1[4] = steps*interdata[4];
            k1[5] = steps*interdata[5];
            k1[6] = steps*interdata[6];
            k1[7] = steps*interdata[7];

            interdata[0] = interdata[0]+0.5*k1[4];
            interdata[1] = interdata[1]+0.5*k1[5];
            interdata[2] = interdata[2]+0.5*k1[6];
            interdata[3] = interdata[3]+0.5*k1[7];
            interdata[4] = interdata[4]+0.5*k1[0];
            interdata[5] = interdata[5]+0.5*k1[1];
            interdata[6] = interdata[6]+0.5*k1[2];
            interdata[7] = interdata[7]+0.5*k1[3];
            interdata[8] = G*m2*/(((((interdata[0]-interdata[2])^2)+((interdata[1]-interdata[3])^2))^(1/2))^3)*(interdata[2]-interdata[0]);
            interdata[9] = G*m2*/(((((interdata[0]-interdata[2])^2)+((interdata[1]-interdata[3])^2))^(1/2))^3)*(interdata[3]-interdata[1]);
            interdata[10] = G*m1*/(((((interdata[0]-interdata[2])^2)+((interdata[1]-interdata[3])^2))^(1/2))^3)*(interdata[0]-interdata[2]);
            interdata[11] = G*m1*/(((((interdata[0]-interdata[2])^2)+((interdata[1]-interdata[3])^2))^(1/2))^3)*(interdata[1]-interdata[3]);

            //k2
           
            k2[0] = steps*interdata[8];
            k2[1] = steps*interdata[9];
            k2[2] = steps*interdata[10];
            k2[3] = steps*interdata[11];
            k2[4] = steps*interdata[4];
            k2[5] = steps*interdata[5];
            k2[6] = steps*interdata[6];
            k2[7] = steps*interdata[7];

            interdata[0] = interdata[0]+0.5*k2[4];
            interdata[1] = interdata[1]+0.5*k2[5];
            interdata[2] = interdata[2]+0.5*k2[6];
            interdata[3] = interdata[3]+0.5*k2[7];
            interdata[4] = interdata[4]+0.5*k2[0];
            interdata[5] = interdata[5]+0.5*k2[1];
            interdata[6] = interdata[6]+0.5*k2[2];
            interdata[7] = interdata[7]+0.5*k2[3];
            interdata[8] = G*m2*/(((((interdata[0]-interdata[2])^2)+((interdata[1]-interdata[3])^2))^(1/2))^3)*(interdata[2]-interdata[0]);
            interdata[9] = G*m2*/(((((interdata[0]-interdata[2])^2)+((interdata[1]-interdata[3])^2))^(1/2))^3)*(interdata[3]-interdata[1]);
            interdata[10] = G*m1*/(((((interdata[0]-interdata[2])^2)+((interdata[1]-interdata[3])^2))^(1/2))^3)*(interdata[0]-interdata[2]);
            interdata[11] = G*m1*/(((((interdata[0]-interdata[2])^2)+((interdata[1]-interdata[3])^2))^(1/2))^3)*(interdata[1]-interdata[3]);

            //k3
           
            k3[0] = steps*interdata[8];
            k3[1] = steps*interdata[9];
            k3[2] = steps*interdata[10];
            k3[3] = steps*interdata[11];
            k3[4] = steps*interdata[4];
            k3[5] = steps*interdata[5];
            k3[6] = steps*interdata[6];
            k3[7] = steps*interdata[7];

            interdata[0] = interdata[0]+k3[4];
            interdata[1] = interdata[1]+k3[5];
            interdata[2] = interdata[2]+k3[6];
            interdata[3] = interdata[3]+k3[7];
            interdata[4] = interdata[4]+k3[0];
            interdata[5] = interdata[5]+k3[1];
            interdata[6] = interdata[6]+k3[2];
            interdata[7] = interdata[7]+k3[3];
            interdata[8] = G*m2*/(((((interdata[0]-interdata[2])^2)+((interdata[1]-interdata[3])^2))^(1/2))^3)*(interdata[2]-interdata[0]);
            interdata[9] = G*m2*/(((((interdata[0]-interdata[2])^2)+((interdata[1]-interdata[3])^2))^(1/2))^3)*(interdata[3]-interdata[1]);
            interdata[10] = G*m1*/(((((interdata[0]-interdata[2])^2)+((interdata[1]-interdata[3])^2))^(1/2))^3)*(interdata[0]-interdata[2]);
            interdata[11] = G*m1*/(((((interdata[0]-interdata[2])^2)+((interdata[1]-interdata[3])^2))^(1/2))^3)*(interdata[1]-interdata[3]);

            //k4
           
            k4[0] = steps*interdata[8];
            k4[1] = steps*interdata[9];
            k4[2] = steps*interdata[10];
            k4[3] = steps*interdata[11];
            k4[4] = steps*interdata[4];
            k4[5] = steps*interdata[5];
            k4[6] = steps*interdata[6];
            k4[7] = steps*interdata[7];

            data[0][i+1] = data[0][i]+(1.0/6.0)*(k1[0]+2*k2[0]+2*k3[0]+k4[0]);
            data[1][i+1] = data[1][i]+(1.0/6.0)*(k1[1]+2*k2[1]+2*k3[1]+k4[1]);
            data[2][i+1] = data[2][i]+(1.0/6.0)*(k1[2]+2*k2[2]+2*k3[2]+k4[2]);
            data[3][i+1] = data[3][i]+(1.0/6.0)*(k1[3]+2*k2[3]+2*k3[3]+k4[3]);
            data[4][i+1] = data[4][i]+(1.0/6.0)*(k1[4]+2*k2[4]+2*k3[4]+k4[4]);
            data[5][i+1] = data[5][i]+(1.0/6.0)*(k1[5]+2*k2[5]+2*k3[5]+k4[5]);
            data[6][i+1] = data[6][i]+(1.0/6.0)*(k1[6]+2*k2[6]+2*k3[6]+k4[6]);
            data[7][i+1] = data[7][i]+(1.0/6.0)*(k1[7]+2*k2[7]+2*k3[7]+k4[7]);

        }
        //output new values
    for (int i=0; i<datalength; i++){
    outfile1<<i*steps<<'\t'<<data[0][i]<<'\t'<<data[1][i]<<'\t'<<data[4][i]<<'\t'<<data[5][i]<<'\t'<<data[2][i]<<'\t'<<data[3][i]<<'\t'<<data[6][i]<<'\t'<<data[7][i]<<endl;
    }
        outfile1.close();
    }
     
     
  2. jcsd
  3. Dec 14, 2012 #2

    SteamKing

    User Avatar
    Staff Emeritus
    Science Advisor
    Homework Helper

    I'm no C++ expert, but I think there's something funny when void RK4_func has a variable list which consists of 'double,double,double, double, etc.)

    Shouldn't there be some actual parameter names in there somewhere?
     
  4. Dec 14, 2012 #3
    Forward declarations don't need them. The actual function is much farther down.
     
  5. Dec 14, 2012 #4
    OP, aren't you getting errors in your compiler? What do they say? You should be, because I'm getting a ton.
     
  6. Dec 14, 2012 #5

    jtbell

    User Avatar

    Staff: Mentor

    Don't you think it would help if you tell us the error messages that the compiler is giving you? :confused:
     
  7. Dec 15, 2012 #6
    You have lots of typos and some more principal problems.
    Always start from the top of your compiler error output. Correct simple typos (commas, semicolons, etc.)
    Exponentiation is different in c++, you need to import pow() function or overload ^ as such beforehand.
     
  8. Dec 15, 2012 #7
    Sorry sorry, yes I've made some fundamental problems. Not very used to C++, I'm rewriting it now and I'll post up the new code. I'll also be deleting the forward declaration, and making it easier to read.
     
  9. Dec 15, 2012 #8

    AlephZero

    User Avatar
    Science Advisor
    Homework Helper

    People often put variable names in function declarations, but strictly speaking it's better not to.

    Think what could happen if a "variable name" was redefined in a header file that you included, for example. At best you would get a very confusing error message. At worst, the code would be valid, but not what you thought you had written.
     
  10. Dec 15, 2012 #9
    I redid it, but I still get this error.

    TwoBodies.cpp: In function ‘void RK4_func(double, double, double, double, double, double, double, double, double, double, double, double)’:
    TwoBodies.cpp:87: error: size of array ‘data’ has non-integral type ‘double’

    line 87 is double data [8][datalength];

    The code is
    Code (Text):
    #include <iostream>
    #include <vector>
    #include <fstream>
    #include <math.h>

    using namespace std;

    void RK4_func (double, double, double, double, double, double, double, double, double, double, double, double);

    int main(){

        //initialise variables
        double time_end,time_step,mass1,x1,y1,vx1,vy1,mass2,x2,y2,vx2,vy2;

        //Get the value into variables
        cout << "Enter Simulation Duration in seconds" << endl;
        cin >> time_end;
        cout << "\n";

        cout << "Enter Precision (Time step. Recommended 100 to 300 seconds )" << endl;
        cin >> time_step;
        cout << "\n";

        cout << "Enter mass of body 1 with 10^20kg scale" << endl;
        cin >> mass1;
        cout << "\n";

        cout << "Enter initial body 1 x position" << endl;
        cin >> x1;
        cout << "\n";

        cout << "Enter initial body 1 y position" << endl;
        cin >> y1;
        cout << "\n";

        cout << "Enter initial body 1 x velocity" << endl;
        cin >> vx1;
        cout << "\n";

        cout << "Enter initial body 1 x velocity" << endl;
        cin >> vy1;
        cout << "\n";

        cout << "Enter mass of body 2 with 10^20kg scale" << endl;
        cin >> mass2;
        cout << "\n";

        cout << "Enter initial body 2 x position" << endl;
        cin >> x2;
        cout << "\n";

        cout << "Enter initial body 2 y position" << endl;
        cin >> y2;
        cout << "\n";

        cout << "Enter initial body 2 x velocity" << endl;
        cin >> vx2;
        cout << "\n";

        cout << "Enter initial body 2 y velocity" << endl;
        cin >> vy2;
        cout << "\n";


        //Pass values into RK4_func and Output data via Excel File

        RK4_func(time_end,time_step,mass1,x1,y1,vx1,vy1,mass2,x2,y2,vx2,vy2);

        return 0;
    }

    //Assume function enters data in SI units

    void RK4_func (double timeend, double steps, double m1, double x1, double y1, double vx1, double vy1, double m2, double x2, double y2, double vx2,double vy2)
    {

        //Excel file creation

        ofstream outfile1 ("RK4.xls");
        outfile1<<"time"<<'\t'<<"x1"<<'\t'<<"y1"<<'\t'<<"vx1"<<'\t'<<"vy1"<<'\t'<<"x2"<<'\t'<<"y2"<<'\t'<<"vx2"<<'\t'<<"vy2"<<endl;
        outfile1<<0<<'\t'<<x1<<'\t'<<y1<<'\t'<<vx1<<'\t'<<vy1<<'\t'<<x2<<'\t'<<y2<<'\t'<<vx2<<'\t'<<vy2<<endl;

        //Declare constants and array

        const double G = 6.67*pow(10,-11);
        double datalength = 1+(timeend/steps);
        double data[8][datalength];

        data[0][0]=x1;
        data[1][0]=y1;
        data[2][0]=x2;
        data[3][0]=y2;
        data[4][0]=vx1;
        data[5][0]=vy1;
        data[6][0]=vx2;
        data[7][0]=vy2;

        for (int i=0; i<datalength; i++)
        {
            vector<double> interdata(12), k1 (8), k2 (8), k3 (8), k4 (8);

            // let 0~3 indicate position x1y1x2y2 4~7 indicate vx1vy1vx2vy2 8~11 indicate ax1ay1ax2ay2

            //acceleration x1 y1 x2 y2

            interdata[0]=data[0][i];
            interdata[1]=data[1][i];
            interdata[2]=data[2][i];
            interdata[3]=data[3][i];
            interdata[4]=data[4][i];
            interdata[5]=data[5][i];
            interdata[6]=data[6][i];
            interdata[7]=data[7][i];
            interdata[8] = G*m2/pow(pow((pow((data[0][i]-data[2][i]),2)+pow((data[1][i]-data[3][i]),2)),(1/2)),3)*(data[2][i]-data[0][i]);
            interdata[9] = G*m2/pow(pow((pow((data[0][i]-data[2][i]),2)+pow((data[1][i]-data[3][i]),2)),(1/2)),3)*(data[3][i]-data[1][i]);
            interdata[10] = G*m1/pow(pow((pow((data[0][i]-data[2][i]),2)+pow((data[1][i]-data[3][i]),2)),(1/2)),3)*(data[0][i]-data[2][i]);
            interdata[11] = G*m1/pow(pow((pow((data[0][i]-data[2][i]),2)+pow((data[1][i]-data[3][i]),2)),(1/2)),3)*(data[1][i]-data[3][i]);

            //k1 is Euler Step

            k1[0] = steps*interdata[8];
            k1[1] = steps*interdata[9];
            k1[2] = steps*interdata[10];
            k1[3] = steps*interdata[11];
            k1[4] = steps*interdata[4];
            k1[5] = steps*interdata[5];
            k1[6] = steps*interdata[6];
            k1[7] = steps*interdata[7];

            interdata[0] = interdata[0]+0.5*k1[4];
            interdata[1] = interdata[1]+0.5*k1[5];
            interdata[2] = interdata[2]+0.5*k1[6];
            interdata[3] = interdata[3]+0.5*k1[7];
            interdata[4] = interdata[4]+0.5*k1[0];
            interdata[5] = interdata[5]+0.5*k1[1];
            interdata[6] = interdata[6]+0.5*k1[2];
            interdata[7] = interdata[7]+0.5*k1[3];
            interdata[8] = G*m2/pow(pow((pow((interdata[0]-interdata[2]),2)+pow((interdata[1]-interdata[3]),2)),(1/2)),3)*(interdata[2]-interdata[0]);
            interdata[9] = G*m2/pow(pow((pow((interdata[0]-interdata[2]),2)+pow((interdata[1]-interdata[3]),2)),(1/2)),3)*(interdata[3]-interdata[1]);
            interdata[10] = G*m1/pow(pow((pow((interdata[0]-interdata[2]),2)+pow((interdata[1]-interdata[3]),2)),(1/2)),3)*(interdata[0]-interdata[2]);
            interdata[11] = G*m1/pow(pow((pow((interdata[0]-interdata[2]),2)+pow((interdata[1]-interdata[3]),2)),(1/2)),3)*(interdata[1]-interdata[3]);

            //k2

            k2[0] = steps*interdata[8];
            k2[1] = steps*interdata[9];
            k2[2] = steps*interdata[10];
            k2[3] = steps*interdata[11];
            k2[4] = steps*interdata[4];
            k2[5] = steps*interdata[5];
            k2[6] = steps*interdata[6];
            k2[7] = steps*interdata[7];

            interdata[0] = interdata[0]+0.5*k2[4];
            interdata[1] = interdata[1]+0.5*k2[5];
            interdata[2] = interdata[2]+0.5*k2[6];
            interdata[3] = interdata[3]+0.5*k2[7];
            interdata[4] = interdata[4]+0.5*k2[0];
            interdata[5] = interdata[5]+0.5*k2[1];
            interdata[6] = interdata[6]+0.5*k2[2];
            interdata[7] = interdata[7]+0.5*k2[3];
            interdata[8] = G*m2/pow(pow((pow((interdata[0]-interdata[2]),2)+pow((interdata[1]-interdata[3]),2)),(1/2)),3)*(interdata[2]-interdata[0]);
            interdata[9] = G*m2/pow(pow((pow((interdata[0]-interdata[2]),2)+pow((interdata[1]-interdata[3]),2)),(1/2)),3)*(interdata[3]-interdata[1]);
            interdata[10] = G*m1/pow(pow((pow((interdata[0]-interdata[2]),2)+pow((interdata[1]-interdata[3]),2)),(1/2)),3)*(interdata[0]-interdata[2]);
            interdata[11] = G*m1/pow(pow((pow((interdata[0]-interdata[2]),2)+pow((interdata[1]-interdata[3]),2)),(1/2)),3)*(interdata[1]-interdata[3]);

            //k3

            k3[0] = steps*interdata[8];
            k3[1] = steps*interdata[9];
            k3[2] = steps*interdata[10];
            k3[3] = steps*interdata[11];
            k3[4] = steps*interdata[4];
            k3[5] = steps*interdata[5];
            k3[6] = steps*interdata[6];
            k3[7] = steps*interdata[7];

            interdata[0] = interdata[0]+k3[4];
            interdata[1] = interdata[1]+k3[5];
            interdata[2] = interdata[2]+k3[6];
            interdata[3] = interdata[3]+k3[7];
            interdata[4] = interdata[4]+k3[0];
            interdata[5] = interdata[5]+k3[1];
            interdata[6] = interdata[6]+k3[2];
            interdata[7] = interdata[7]+k3[3];
            interdata[8] = G*m2/pow(pow((pow((interdata[0]-interdata[2]),2)+pow((interdata[1]-interdata[3]),2)),(1/2)),3)*(interdata[2]-interdata[0]);
            interdata[9] = G*m2/pow(pow((pow((interdata[0]-interdata[2]),2)+pow((interdata[1]-interdata[3]),2)),(1/2)),3)*(interdata[3]-interdata[1]);
            interdata[10] = G*m1/pow(pow((pow((interdata[0]-interdata[2]),2)+pow((interdata[1]-interdata[3]),2)),(1/2)),3)*(interdata[0]-interdata[2]);
            interdata[11] = G*m1/pow(pow((pow((interdata[0]-interdata[2]),2)+pow((interdata[1]-interdata[3]),2)),(1/2)),3)*(interdata[1]-interdata[3]);

            //k4

            k4[0] = steps*interdata[8];
            k4[1] = steps*interdata[9];
            k4[2] = steps*interdata[10];
            k4[3] = steps*interdata[11];
            k4[4] = steps*interdata[4];
            k4[5] = steps*interdata[5];
            k4[6] = steps*interdata[6];
            k4[7] = steps*interdata[7];

            //Form the next position and velocity
           
            data[0][i+1] = data[0][i]+(1.0/6.0)*(k1[0]+2*k2[0]+2*k3[0]+k4[0]);
            data[1][i+1] = data[1][i]+(1.0/6.0)*(k1[1]+2*k2[1]+2*k3[1]+k4[1]);
            data[2][i+1] = data[2][i]+(1.0/6.0)*(k1[2]+2*k2[2]+2*k3[2]+k4[2]);
            data[3][i+1] = data[3][i]+(1.0/6.0)*(k1[3]+2*k2[3]+2*k3[3]+k4[3]);
            data[4][i+1] = data[4][i]+(1.0/6.0)*(k1[4]+2*k2[4]+2*k3[4]+k4[4]);
            data[5][i+1] = data[5][i]+(1.0/6.0)*(k1[5]+2*k2[5]+2*k3[5]+k4[5]);
            data[6][i+1] = data[6][i]+(1.0/6.0)*(k1[6]+2*k2[6]+2*k3[6]+k4[6]);
            data[7][i+1] = data[7][i]+(1.0/6.0)*(k1[7]+2*k2[7]+2*k3[7]+k4[7]);

        }
        //output new values
        for (int i=0; i<datalength; i++){
            outfile1<<i*steps<<'\t'<<data[0][i]<<'\t'<<data[1][i]<<'\t'<<data[4][i]<<'\t'<<data[5][i]<<'\t'<<data[2][i]<<'\t'<<data[3][i]<<'\t'<<data[6][i]<<'\t'<<data[7][i]<<endl;
        }
        outfile1.close();
    }

    // Display via GnuPlot
     
     
    Last edited: Dec 15, 2012
  11. Dec 15, 2012 #10
    You must used a fixed point type as an index into an array, you can't use a floating point value. Likewise you must have a fixed point count of objects (same thing). You can't have 3.2 elements in the array...
     
  12. Dec 15, 2012 #11
    so how do i specify that? do i use long int or int?

    EDIT: Yes. I switched the data into unsigned long int, then it compiled! but now the RK4 seems off since it doesn't give me a result.
     
  13. Dec 15, 2012 #12
    You have no experience in programming? It sounds like you don't even know where to begin when it comes to debugging.
     
  14. Dec 15, 2012 #13
    Sorry, I don't really have that much experience. It's my first course, and it's only been a few months since I started learning C++/programming. But it's compiling now. I just think the RK4 is off now, because it gives me the excel file, but the numbers are ridiculous. Is there something that I did wrong with my RK4? I mean, i'm sure there must be, but I don't see it.

    The RK4 function is here:
    Code (Text):
    void RK4_func (double timeend, double steps, double m1, double x1, double y1, double vx1, double vy1, double m2, double x2, double y2, double vx2,double vy2)
    {

        //Excel file creation

        ofstream outfile1 ("RK4.xls");
        outfile1<<"time"<<'\t'<<"x1"<<'\t'<<"y1"<<'\t'<<"vx1"<<'\t'<<"vy1"<<'\t'<<"x2"<<'\t'<<"y2"<<'\t'<<"vx2"<<'\t'<<"vy2"<<endl;
        outfile1<<0<<'\t'<<x1<<'\t'<<y1<<'\t'<<vx1<<'\t'<<vy1<<'\t'<<x2<<'\t'<<y2<<'\t'<<vx2<<'\t'<<vy2<<endl;

        //Declare constants and array

        const double G = 6.67*pow(10,-11);
        unsigned long int datalength = 1+(timeend/steps);
        double data[8][datalength];

        data[0][0]=x1;
        data[1][0]=y1;
        data[2][0]=x2;
        data[3][0]=y2;
        data[4][0]=vx1;
        data[5][0]=vy1;
        data[6][0]=vx2;
        data[7][0]=vy2;

        for (int i=0; i<datalength; i++)
        {
            vector<double> interdata(12), k1 (8), k2 (8), k3 (8), k4 (8);

            // let 0~3 indicate position x1y1x2y2 4~7 indicate vx1vy1vx2vy2 8~11 indicate ax1ay1ax2ay2

            //acceleration x1 y1 x2 y2

            interdata[0]=data[0][i];
            interdata[1]=data[1][i];
            interdata[2]=data[2][i];
            interdata[3]=data[3][i];
            interdata[4]=data[4][i];
            interdata[5]=data[5][i];
            interdata[6]=data[6][i];
            interdata[7]=data[7][i];
            interdata[8] = G*m2/pow(pow((pow((data[0][i]-data[2][i]),2)+pow((data[1][i]-data[3][i]),2)),(1/2)),3)*(data[2][i]-data[0][i]);
            interdata[9] = G*m2/pow(pow((pow((data[0][i]-data[2][i]),2)+pow((data[1][i]-data[3][i]),2)),(1/2)),3)*(data[3][i]-data[1][i]);
            interdata[10] = G*m1/pow(pow((pow((data[0][i]-data[2][i]),2)+pow((data[1][i]-data[3][i]),2)),(1/2)),3)*(data[0][i]-data[2][i]);
            interdata[11] = G*m1/pow(pow((pow((data[0][i]-data[2][i]),2)+pow((data[1][i]-data[3][i]),2)),(1/2)),3)*(data[1][i]-data[3][i]);

            //k1 is Euler Step

            k1[0] = steps*interdata[8];
            k1[1] = steps*interdata[9];
            k1[2] = steps*interdata[10];
            k1[3] = steps*interdata[11];
            k1[4] = steps*interdata[4];
            k1[5] = steps*interdata[5];
            k1[6] = steps*interdata[6];
            k1[7] = steps*interdata[7];

            interdata[0] = interdata[0]+0.5*k1[4];
            interdata[1] = interdata[1]+0.5*k1[5];
            interdata[2] = interdata[2]+0.5*k1[6];
            interdata[3] = interdata[3]+0.5*k1[7];
            interdata[4] = interdata[4]+0.5*k1[0];
            interdata[5] = interdata[5]+0.5*k1[1];
            interdata[6] = interdata[6]+0.5*k1[2];
            interdata[7] = interdata[7]+0.5*k1[3];
            interdata[8] = G*m2/pow(pow((pow((interdata[0]-interdata[2]),2)+pow((interdata[1]-interdata[3]),2)),(1/2)),3)*(interdata[2]-interdata[0]);
            interdata[9] = G*m2/pow(pow((pow((interdata[0]-interdata[2]),2)+pow((interdata[1]-interdata[3]),2)),(1/2)),3)*(interdata[3]-interdata[1]);
            interdata[10] = G*m1/pow(pow((pow((interdata[0]-interdata[2]),2)+pow((interdata[1]-interdata[3]),2)),(1/2)),3)*(interdata[0]-interdata[2]);
            interdata[11] = G*m1/pow(pow((pow((interdata[0]-interdata[2]),2)+pow((interdata[1]-interdata[3]),2)),(1/2)),3)*(interdata[1]-interdata[3]);

            //k2

            k2[0] = steps*interdata[8];
            k2[1] = steps*interdata[9];
            k2[2] = steps*interdata[10];
            k2[3] = steps*interdata[11];
            k2[4] = steps*interdata[4];
            k2[5] = steps*interdata[5];
            k2[6] = steps*interdata[6];
            k2[7] = steps*interdata[7];

            interdata[0] = interdata[0]+0.5*k2[4];
            interdata[1] = interdata[1]+0.5*k2[5];
            interdata[2] = interdata[2]+0.5*k2[6];
            interdata[3] = interdata[3]+0.5*k2[7];
            interdata[4] = interdata[4]+0.5*k2[0];
            interdata[5] = interdata[5]+0.5*k2[1];
            interdata[6] = interdata[6]+0.5*k2[2];
            interdata[7] = interdata[7]+0.5*k2[3];
            interdata[8] = G*m2/pow(pow((pow((interdata[0]-interdata[2]),2)+pow((interdata[1]-interdata[3]),2)),(1/2)),3)*(interdata[2]-interdata[0]);
            interdata[9] = G*m2/pow(pow((pow((interdata[0]-interdata[2]),2)+pow((interdata[1]-interdata[3]),2)),(1/2)),3)*(interdata[3]-interdata[1]);
            interdata[10] = G*m1/pow(pow((pow((interdata[0]-interdata[2]),2)+pow((interdata[1]-interdata[3]),2)),(1/2)),3)*(interdata[0]-interdata[2]);
            interdata[11] = G*m1/pow(pow((pow((interdata[0]-interdata[2]),2)+pow((interdata[1]-interdata[3]),2)),(1/2)),3)*(interdata[1]-interdata[3]);

            //k3

            k3[0] = steps*interdata[8];
            k3[1] = steps*interdata[9];
            k3[2] = steps*interdata[10];
            k3[3] = steps*interdata[11];
            k3[4] = steps*interdata[4];
            k3[5] = steps*interdata[5];
            k3[6] = steps*interdata[6];
            k3[7] = steps*interdata[7];

            interdata[0] = interdata[0]+k3[4];
            interdata[1] = interdata[1]+k3[5];
            interdata[2] = interdata[2]+k3[6];
            interdata[3] = interdata[3]+k3[7];
            interdata[4] = interdata[4]+k3[0];
            interdata[5] = interdata[5]+k3[1];
            interdata[6] = interdata[6]+k3[2];
            interdata[7] = interdata[7]+k3[3];
            interdata[8] = G*m2/pow(pow((pow((interdata[0]-interdata[2]),2)+pow((interdata[1]-interdata[3]),2)),(1/2)),3)*(interdata[2]-interdata[0]);
            interdata[9] = G*m2/pow(pow((pow((interdata[0]-interdata[2]),2)+pow((interdata[1]-interdata[3]),2)),(1/2)),3)*(interdata[3]-interdata[1]);
            interdata[10] = G*m1/pow(pow((pow((interdata[0]-interdata[2]),2)+pow((interdata[1]-interdata[3]),2)),(1/2)),3)*(interdata[0]-interdata[2]);
            interdata[11] = G*m1/pow(pow((pow((interdata[0]-interdata[2]),2)+pow((interdata[1]-interdata[3]),2)),(1/2)),3)*(interdata[1]-interdata[3]);

            //k4

            k4[0] = steps*interdata[8];
            k4[1] = steps*interdata[9];
            k4[2] = steps*interdata[10];
            k4[3] = steps*interdata[11];
            k4[4] = steps*interdata[4];
            k4[5] = steps*interdata[5];
            k4[6] = steps*interdata[6];
            k4[7] = steps*interdata[7];

            //Form the next position and velocity

            data[0][i+1] = data[0][i]+(1.0/6.0)*(k1[0]+2*k2[0]+2*k3[0]+k4[0]);
            data[1][i+1] = data[1][i]+(1.0/6.0)*(k1[1]+2*k2[1]+2*k3[1]+k4[1]);
            data[2][i+1] = data[2][i]+(1.0/6.0)*(k1[2]+2*k2[2]+2*k3[2]+k4[2]);
            data[3][i+1] = data[3][i]+(1.0/6.0)*(k1[3]+2*k2[3]+2*k3[3]+k4[3]);
            data[4][i+1] = data[4][i]+(1.0/6.0)*(k1[4]+2*k2[4]+2*k3[4]+k4[4]);
            data[5][i+1] = data[5][i]+(1.0/6.0)*(k1[5]+2*k2[5]+2*k3[5]+k4[5]);
            data[6][i+1] = data[6][i]+(1.0/6.0)*(k1[6]+2*k2[6]+2*k3[6]+k4[6]);
            data[7][i+1] = data[7][i]+(1.0/6.0)*(k1[7]+2*k2[7]+2*k3[7]+k4[7]);

        }
        //output new values
        for (int i=0; i<datalength; i++){
            outfile1<<i*steps<<'\t'<<data[0][i]<<'\t'<<data[1][i]<<'\t'<<data[4][i]<<'\t'<<data[5][i]<<'\t'<<data[2][i]<<'\t'<<data[3][i]<<'\t'<<data[6][i]<<'\t'<<data[7][i]<<endl;
        }
        outfile1.close();
    }
     
    Last edited: Dec 15, 2012
  15. Dec 19, 2012 #14
    Being confused is normal. Knowing how to dissect the problems are the key.

    To be honest, I don't understand your code. k1[0] to k1[3] should give the derivative of coordinates (velocity) times the stepsize, but you have second derivative instead (the acceleration).
    Same goes for k1[4] to k[7]. Shouldn't they have stepsize times the acceleration instead?
    Of course then you seem to flip it back to correct places in interdata, but the final calculation uses k vectors where the problem remains.

    Maybe it's possible to correct the problem by only changing the coefficient vector numbering in
    data[0][i+1] = data[0] + ...
    ...
    ...
    etc

    If all else fails try doing only Euler step method at first ( comment the rest out). Then there's less to think about.
     
Know someone interested in this topic? Share this thread via Reddit, Google+, Twitter, or Facebook




Similar Discussions: Two Body problem with RK4 in C++
  1. Problems in C++ (Replies: 5)

  2. A problem with C (Replies: 4)

Loading...