Problems with Runge-Kutta method

[Leonardo Machado]

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.

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];
          
      
      }
}

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.

 

[D H]

@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.

 

[Leonardo Machado]

@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.

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.

 

[D H]

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

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.

 

[rcgldr]

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

 

[Leonardo Machado]

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.