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

Tags:
1. Jul 24, 2016

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

}
}

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. Jul 25, 2016

### D H

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

3. Jul 25, 2016

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.

4. Jul 25, 2016

### D H

Staff Emeritus
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.

5. Jul 25, 2016

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

6. Jul 26, 2016