Runge-Kutta Method for Second Order Differential Equations

Description

The classical Runge-Kutta method applied to the second order differential equation y''(x) = f(x, y, y') with initial conditions y(x₀) = y₀ and y'(x₀) = y'₀ evaluates the function f(x,y,y') four times per step and can be derived by transforming the problem to a coupled system of first order differential equations.

The algorithm for the implicit central difference method is derived by using the central difference approximations for y'(x) and y''(x):

y'(x) = (y(x+h) - y(x-h) ) / 2h, and
y''(x) = (y(x+h) - 2y(x) + y(x-h)) / h².

Let x_n = x₀ + nh, y_n be the approximation to y(x_n), and f_n = f(x_n, y_n, (y_{n + 1} - y_{n - 1} ) / 2h ), the procedure proceeds recursively via the implicit equation for y_{n + 1} as follows:

y_{n + 1} = 2 y_n - y_{n - 1} + h² f_n.

Function List

Functions:

void Runge_Kutta_2nd_Order( double (*f)(double,double,double), double x0, double y[ ], double c, double h, int number_of_steps )

This function uses the classical Runge-Kutta method to estimate the solution of the initial value problem, y'' = f(x,y,y'); y(x0) = y[0] and y'(x0) = c, at x0 + nh where for n =1, . . ., number_of_steps and h is the step size. On input, y[0] is the value of y(x) at x = x0. On output y[n] is the value of y(x) at x = x0 + n h for n =0, . . ., number_of_steps.

C Source

File: runge_kutta_2nd_order.c

////////////////////////////////////////////////////////////////////////////////
// File: runge_kutta_2nd_order.c                                              //
// Routines:                                                                  //
//    Runge_Kutta_2nd_Order                                                   //
////////////////////////////////////////////////////////////////////////////////

////////////////////////////////////////////////////////////////////////////////
//                                                                            //
//  Description:                                                              //
//     A differential equation of order greater than one can be transformed   //
//     to a multidimensional first order differential equation.  In           //
//     particular, given a second order differential equation y'' = f(x,y,y') //
//     define z = y', then the pair (y,z)' = (z, f(x,y,z) ) is a first order  //
//     differential equation.  The initial conditions y(x0) = y[0],           //
//     y'(x0) = c becomes the initial condition (y,z)(x0) = (y[0],c).         //
//     when x = x0 numerically evaluates f(x,y) four times per step. Then for //
//     step i+1,                                                              //
//           y[i+1] = y[i] + h * ( yp[i] + 1/6 (k1 + k2 + k3 ) and            //
//           yp[i+1] = yp[i] + 1/6 * ( k1 + 2 k2 + 2 k3 + k4 ) where          //
//         k1 = h * f(x[i], y[i], yp[i]),                                     //
//         k2 = h * f(x[i]+h/2, y[i]+(h/2)*yp[i], yp[i]+k1/2),                //
//         k3 = h * f(x[i]+h/2, y[i]+(h/2)*yp[i]+(h/4)*k1, yp[i]+k2/2),       //
//     x[i] = x0 + i * h.                                                     //
//                                                                            //
////////////////////////////////////////////////////////////////////////////////

static const double one_sixth = 1.0 / 6.0;

////////////////////////////////////////////////////////////////////////////////
//  void Runge_Kutta_2nd_Order( double (*f)(double, double, double),          //
//          double x0, double y[], double c, double h, int number_of_steps ); //
//                                                                            //
//  Description:                                                              //
//     This function solves the 2nd order differential equation y''=f(x,y,y') //
//     with the initial conditions y=y[0] and y'= c at x = x0.  The values    //
//     returned are y[n] which is the value of y evaluated at x = x0 + n*h.   //
//                                                                            //
//  Arguments:                                                                //
//     double *f  Pointer to the function which returns the rate of change of //
//                the slope with respect to x at the point (x,y) with slope   //
//                y'(x) of integral curve of the second order differential    //
//                equation y'' = f(x,y,y') which passes through the point     //
//                (x0,y[0]) with initial slope y'(x0) = c.                    //
//     double y[] On input y[0] is the initial value of y at x = x0, on output//
//                for i >= 1,  y[i] is given by the algorithm described above.//
//     double x0  Initial value of x.                                         //
//     double h   Step size                                                   //
//     int    number_of_steps  The number of steps, y[] should be dimensioned //
//                             number_of_steps + 1.                           //
//                                                                            //
//  Return Values:                                                            //
//     This routine is of type void and hence does not return a value.  The   //
//     solution of y'' = f(x,y,y') from x = x0 to x = x0 + number_of_steps*h  //
//     is stored in the input array y[].                                      //
//                                                                            //
////////////////////////////////////////////////////////////////////////////////
//                                                                            //
void Runge_Kutta_2nd_Order( double (*f)(double, double, double), double x0,
                        double y[], double c, double h, int number_of_steps ) {

   double k1, k2, k3, k4;
   double h2 = 0.5 * h;
   double ych2;
   int i;

   for (i = 0; i < number_of_steps; x0 += h, i++) {
      ych2 = y[i] + c * h2;
      k1 = h * (*f)(x0, y[i], c);
      k2 = h * (*f)(x0+h2, ych2, c + 0.5 * k1);
      k3 = h * (*f)(x0+h2, ych2 + 0.25 * k1 * h, c + 0.5 * k2);
      k4 = h * (*f)(x0+h, y[i] + c * h + h2 * k2 , c + k3);
      y[i+1] = y[i] + ( c + one_sixth * (k1 + k2 + k3) ) * h;
      c += one_sixth * ( k1 + k2 + k2 + k3 + k3 + k4 );
   }
}