The Runge Kutta 4th order Algorithm

I am doing an intern at Zeus Numerix where I have to simulate motion of an artillery missile fired from a launch vehicle.After doing the dynamics part of the analysis we get a series of second order differential equations which are solved by the Runge – Kutta method.

The Runge-Kutta algorithm is the magic formula behind many of the physics simulations . The Runge-Kutta algorithm lets us solve a differential equation numerically (that is, approximately) it is known to be very accurate and well-behaved for a wide range of problems.

Consider the single variable problem x‘ = f (t, x) with initial condition x(0) = x0. Suppose that xn is the value of the variable at time tn. The Runge-Kutta formula takes xn and tn and calculates an approximation for xn+1 at a brief time later, tn+h. It uses a weighted average of approximated values of f (t, x) at several times within the interval (tn, tn+h). The formula is given by xn+1 = xn + h6 (a + 2 b + 2 c + d) where a = f (tn, xn)
b = f (tn + h2, xn + h2 a)
c = f (tn + h2, xn + h2 b)
d = f (tn + h, xn + h c).

To run the simulation, we start with x0 and find x1 using the formula above. Then we plug in x1 to find x2 and so on.

With multiple variables, the Runge-Kutta algorithm looks similar to the above equations, except that the variables become vectors.

Make sure that u keep all the derivatives on one side of the equation and then solve.

A C++ code for it is as follows:

double solve(double timestep,double init,double stepsize){
double ret=init;
double time=0.0;
while(timestep–){

double k1,k2,k3,k4;//fc is the function obtained to the right hand side of the diff. equation
k1=fc(time,ret);
k2=fc(time+(stepsize/2.0),ret+(1.0/2.0 * k1 * stepsize));
k3=fc(time+(stepsize/2.0),ret+(1.0/2.0 * k2 * stepsize));
k4=fc(time+(stepsize/2.0),ret+(k3 * stepsize));
ret+=((k1+2*k2+2*k3+k4)/6.0);

}
return ret;
}

Edit: Thanks to Fabio for the correction  in the code