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

C/++/# Problems with Runge-Kutta method

  1. Jul 24, 2016 #1
    I'm writing a program to compute an ODE solution of the Kepler's problem based on Runge-Kutta 4th order method and today i've past the whole day trying to made it work, but i've failed, maybe you could help me to kill the problem ?
    The solutions is cartesian.

    Code (C):
    int main(){
     
        int n;
     
        double tau, x[4][200], xk[4], f[4][200], GM, PI = 3.14159265;
     
        double f1[4], f2[4], f3[4], f4[4] ;
     
        GM = 4 * pow ( PI, 2);
     
        x[2][1]= 0;
     
        x[3][1]= 0;
     
        cout << "Insert the inicial distance ( AU): " ;
        cin >> x[1][1];
     
        cout << "Choose tau ( Years): " ;
        cin >> tau;
     
        //Circular mov condition
     
        x[4][1]= sqrt ( GM / sqrt ( pow ( x[1][1], 2) + pow ( x[2][1], 2)));
     

        for ( n=1 ; n <= 199 ; n++) {
         
            // Definitions of runge-kutta parametres F1, F2, F3 e F4
         
            xk[1]=0;
            xk[2]=0;
            xk[3]=0;
            xk[4]=0;
         
            f[1][n]= x[3][n] + xk[3];
         
            f[2][n]= x[4][n] + xk[4];
         
            f[3][n]= -GM * (x[1][n] + xk[1]) / pow ( sqrt (  pow ( x[1][n] + xk[1], 2) + pow ( x[2][n] + xk[2], 2)), 3);
         
            f[4][n]= -GM * (x[2][n] + xk[2]) / pow ( sqrt (  pow ( x[1][n] + xk[1], 2) + pow ( x[2][n] + xk[2], 2)), 3);
         
         
            f1[1]= f[1][n];
            f1[2]= f[2][n];
            f1[3]= f[3][n];
            f1[4]= f[4][n];
         
         
            xk[1]=0.5*tau*f1[1];
            xk[2]=0.5*tau*f1[2]/2;
            xk[3]=0.5*f1[3];
            xk[4]=0.5*tau*f1[4];
         
            f2[1]= f[1][n];
            f2[2]= f[2][n];
            f2[3]= f[3][n];
            f2[4]= f[4][n];
         
            xk[1]=0.5*tau*f2[1];
            xk[2]=0.5*tau*f2[2];
            xk[3]=0.5*tau*f2[3];
            xk[4]=0.5*tau*f2[4];
         
            f3[1]= f[1][n];
            f3[2]= f[2][n];
            f3[3]= f[3][n];
            f3[4]= f[4][n];
         
            xk[1]=0.5*tau*f3[1];
            xk[2]=0.5*tau*f3[2];
            xk[3]=0.5*tau*f3[3];
            xk[4]=0.5*tau*f3[4];
               
            f4[1]= f[1][n];
            f4[2]= f[2][n];
            f4[3]= f[3][n];
            f4[4]= f[4][n];
         
            xk[1]=0;
            xk[2]=0;
            xk[3]=0;
            xk[4]=0;
         

            x[1][n+1]= x[1][n] + tau/6*(f1[1]+f4[1]+2*(f2[1]+f3[1]));
            x[2][n+1]= x[2][n] + tau/6*(f1[2]+f4[2]+2*(f2[2]+f3[2]));
            x[3][n+1]= x[3][n] + tau/6*(f1[3]+f4[3]+2*(f2[3]+f3[3]));
            x[4][n+1]= x[4][n] + tau/6*(f1[4]+f4[4]+2*(f2[4]+f3[4]));

            cout << endl
                << "t= " << setw (6) << tau*(n-1)
                << setw (10) << "x= " << setw (6) << x[1][n]
                << setw (10) << "y= " << setw (6) << x[2][n]
                << setw (10) << "vx= " << setw (6) << f[1][n]
                << setw (10) << "vy= " << setw (6) << f[2][n]
                << setw (10) << "fx= " << setw (6) << f[3][n]
                << setw (10) << "fy= " << setw (6) << f[4][n];
             
         
          }
    }
    Erro para o forum.png
    Look like when it takes n=1 there is some error that make the force in x and y equals 0, as you can see in this print screen, also it crashes on n=5 :c

    I would be gratefull if you help me.
     
    Last edited by a moderator: Jul 25, 2016
  2. jcsd
  3. Jul 25, 2016 #2

    D H

    User Avatar
    Staff Emeritus
    Science Advisor

    @Leonardo Machado, I fixed your post to use code tags rather than a bold font.

    Your problem is simple: Arrays in the C family of languages start at zero rather than one.
     
  4. Jul 25, 2016 #3
    Thanks, I'm newbie in here, sorry about the wrong format.



    Dont know if i got it , but when you say it starts at zero you meant for example:

    double x[4] = x[0] , x[1], x[2], x[3] and x[4], right ??

    Case yes, there should not have any problem, cause i'm using n starting at 1 in the loop, these coordinates still exist.
     
  5. Jul 25, 2016 #4

    D H

    User Avatar
    Staff Emeritus
    Science Advisor

    That is not right. What double x[4] does do is to declare x as an array with four elements, accessed as x[0], x[1], x[2], and x[3]. The element x[4] does not exist.

    Your code is incorrect in a number of other ways. You are not implementing RK4 correctly.
     
  6. Jul 25, 2016 #5

    rcgldr

    User Avatar
    Homework Helper

    Example pseudo code RK4 step, p = position, v = velocity, a = acceleration. p[] and v[] don't have to be arrays, just using the syntax of [ i ] as a step indicator. In this case (Kepler), acceleration is a function of position (the velocity parameter is not needed).

    p1 = p[ i ]
    v1 = v[ i ]
    a1 = f(v1, p1)

    p2 = p[ i ] + 1/2 Δt v1
    v2 = v[ i ] + 1/2 Δt a1
    a2 = f(v2, p2)

    p3 = p[ i ] + 1/2 Δt v2
    v3 = v[ i ] + 1/2 Δt a2
    a3 = f(v3, p3)

    p4 = p[ i ] + Δt v3
    v4 = v[ i ] + Δt a3
    a4 = f(v4, p4)

    p[ i+1 ] = p[ i ] + 1/6 Δt (v1 + 2 v2 + 2 v3 + v4)
    v[ i+1 ] = v[ i ] + 1/6 Δt (a1 + 2 a2 + 2 a3 + a4)
    i = i + 1
     
  7. Jul 26, 2016 #6
    I've fixed it, thanks all of you for the contribution and ideas.

    In the end, the mistake was computational, not mathematical, was something about the order that the function was presented.
     
Know someone interested in this topic? Share this thread via Reddit, Google+, Twitter, or Facebook

Have something to add?
Draft saved Draft deleted



Similar Discussions: Problems with Runge-Kutta method
Loading...