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) = *x*_{0}. Suppose that *x _{n}* is the value of the variable at time

*t*. The Runge-Kutta formula takes

_{n}*x*and

_{n}*t*and calculates an approximation for

_{n}*x*

_{n+1}at a brief time later,

*t*+

_{n}*h*. It uses a weighted average of approximated values of

*f*(

*t*,

*x*) at several times within the interval (

*t*,

_{n}*t*+

_{n}*h*). The formula is given by

*x*

_{n+1}=

*x*+

_{n}^{h}⁄

_{6}(

*a*+ 2

*b*+ 2

*c*+

*d*) where

*a*=

*f*(

*t*,

_{n}*x*)

_{n}*b*=

*f*(

*t*+

_{n}^{h}⁄

_{2},

*x*+

_{n}^{h}⁄

_{2}

*a*)

*c*=

*f*(

*t*+

_{n}^{h}⁄

_{2},

*x*+

_{n}^{h}⁄

_{2}

*b*)

*d*=

*f*(

*t*+

_{n}*h*,

*x*+

_{n}*h*

*c*).

To run the simulation, we start with *x*_{0} and find *x*_{1} using the formula above. Then we plug in *x*_{1} to find *x*_{2} 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

nice post jus when i needed it

I think that last evaluating should be:

k4=fc(time+(stepsize),ret+(k3 * stepsize));

по моему мнению: шикарно..

thanks, but i would like to a Fortran code for it.

then do write one …..