Runge-Kutta Integration

Most anybody that has done numerical integration is familiar with Runge Kutta methods. The most common being the fourth order integration equations. But, the equations for simultaneous differential equations are generally not presented so I've put them here.

I'm not trying to teach Runge Kutta methods here, just present the equations for those familiar with numerical integration methods. If this is "greek" to you, take a look at any number of books on numerical methods.

General Runge Kutta Methods

I'm presenting the common second, third and fourth order equations here for completeness. You should be able to find them just about anywhere. These were taken from Handbook of Mathematical Functions (see below). NOTE: there are other equations for the Runge Kutta method, usually differing by constants. I'm not claiming that these are the only methods or even the best.

All equations are for single integration step of the function y = F(x,y) over a step size of h.

Second order
yn+1 = yn + ( k1 + k2 ) / 2 + O( h3 )
k1 = h F( xn,yn )
k2 = h F( xn + h,yn + k1 )
Third order
yn+1 = yn + 1/6 k1 + 2/3 k2 + 1/6 k3 + O( h4 )
k1 = h F( xn,yn )
k2 = h F( xn + h/2, yn + k1/2 )
k3 = h F( xn + h, yn -k1 + 2 k2 )
Fourth order:
yn+1 = yn + 1/6 k1 + 1/3 k2 + 1/3 k3 + 1/6 k4 + O( h5 )
k1 = h F( xn,yn )
k2 = h F( xn + h/2, yn + k1/2 )
k3 = h F( xn + h/2, yn + k2/2 )
k4 = h F( xn + h, yn + k3 )

The equations of motion for a bullet in flight, are simultaneous coupled differential equations. Specifically three equations for acceleration, position and spin. These equations may have been transformed using the substitution:

du/dt = du/dy dy/dt = du/dy uy

producing equations as a function of mach number and range:

a = du/dt = du/dy uy = -p pi/8 CD(m)/SD V V
du/dy = -p pi/8 CD(m)/SD V/uy V

where v = u - w and u is the velocity relative to the ground and w is the wind velocity.

Integrating this equation is much easier when trajectory parameters are required as a function of range. [This vector equation for velocity is derived in the discussion of CD and KD] But, it does introduce another differential equation for time.

Simultaneous Equations of First Order

The Runge Kutta methods for coupled equations are slightly different and not seen very often (at least by me!). I've found them listed in one book, Ballistics of the Future (see below), and have listed them here for anybody that does not have access to this book (only fourth order was shown).

For two equations, dy/dx = F(x, y, z) and dz/dx = G(x, y, z), a single integration step of h is defined by:

yn+1 = yn + 1/6 k1 + 1/3 k2 + 1/3 k3 + 1/6 k4 + O( h5 )
zn+1 = zn + 1/6 l1 + 1/3 l2 + 1/3 l3 + 1/6 l4 + O( h5 )
k1 = h F( xn,yn, zn )
k2 = h F( xn + h/2, yn + k1/2, zn + l1/2 )
k3 = h F( xn + h/2, yn + k2/2, zn + l2/2 )
k4 = h F( xn + h, yn + k3, zn + l3 )
l1 = h G( xn,yn, zn )
l2 = h G( xn + h/2, yn + k1/2, zn + l1/2 )
l3 = h G( xn + h/2, yn + k2/2, zn + l2/2 )
l4 = h G( xn + h, yn + k3, zn + l3 )

Simultaneous Equations of Second Order

Consider the coupled differential equations of second order:

y'' = F ( x, y, z, y', z' )
z'' = G ( x, y, z, y', z' )

where ' denotes differentiation with respect to x.

Assume

y' = u
z' = v

then

u' = F ( x, y, z, u, v )
v' = G ( x, y, z, u, v )

reducing the second order differential equations to a system of four simultaneous eqations of first order.

A single integration step, h of x, is then given by:

yn+1 = yn + 1/6 k1 + 1/3 k2 + 1/3 k3 + 1/6 k4 + O( h5 )
zn+1 = zn + 1/6 l1 + 1/3 l2 + 1/3 l3 + 1/6 l4 + O( h5 )
k1 = h un
k2 = h ( un + m1/2 )
k3 = h ( un + m2/2 )
k4 = h ( un + m3 )
l1 = h vn
l2 = h ( vn + p1/2 )
l3 = h ( vn + p2/2 )
l4 = h ( vn + p3 )
un+1 = un + 1/6 m1 + 1/3 m2 + 1/3 m3 + 1/6 m4 + O( h5 )
vn+1 = vn + 1/6 p1 + 1/3 p2 + 1/3 p3 + 1/6 p4 + O( h5 )
m1 = h F ( xn, yn, zn, un, vn )
m2 = h F ( xn + h/2, yn + k1/2 , zn + l1/2 , un + m1/2 , vn + p1/2 )
m3 = h F ( xn + h/2, yn + k2/2 , zn + l2/2 , un + m2/2 , vn + p2/2 )
m4 = h F ( xn + h, yn + k3 , zn + l3 , un + m3 , vn + p3 )
p1 = h G ( xn, yn, zn, un, vn )
p2 = h G ( xn + h/2, yn + k1/2 , zn + l1/2 , un + m1/2 , vn + p1/2 )
p3 = h G ( xn + h/2, yn + k2/2 , zn + l2/2 , un + m2/2 , vn + p2/2 )
p4 = h G ( xn + h, yn + k3 , zn + l3 , un + m3 , vn + p3 )

References

Ballistic of the Future, J. M. J. Kooy and J. W. H. Uytenbogaart, McGraw-Hill Book Company, Inc. New York, 1946.

Handbook of Mathematical Functions, Edited by Milton Abramowitz and Irene A. Stegun, Dover Publications, Inc., New York, 1972.