Midpoint Method

Description

The midpoint method is an explicit method for approximating the solution of the initial value problem y' = f(x,y); y(x₀) = y₀ at x for a given step size h. For the midpoint method the derivative of y(x) is approximated by the symmetric difference
y'(x) = ( y(x+h) - y(x-h) ) / 2h + O(h²). Then the differential equation becomes
y(x+h) = y(x-h) + 2h f(x,y) + (h³).The approximation y_n for y(x₀+nh) is then given recursively by
y_n+1 = y_n-1 + h f(x_n,y_n)for n = 1, 2, ... .

Locally the midpoint method is a third order method and therefore globally a second order method.

The midpoint method is a stable and convergent method but it is only weakly stable, small perturbations in the initial conditions give rise to growing oscillations.

As a rule, use of the midpoint method should be avoided.

Function List

Functions:

void Midpoint_Method( double (*f)(double, double), double y[ ], double a, double h, int number_of_steps )

This function uses the midpoint method to return the estimate of the solution of the initial value problem, y' = f(x,y); y = y[0] when x = a, at a + nh where for 1 < n < number_of_steps and h is the step size. On input, y[0] is the value of y(x) at x = a and y[1] is the value of y(x) at x = a + h. On output y[n] is the value of y(x) at x = a + n h for 1 < n < number_of_steps.

C Source

File: midpoint_method.c

////////////////////////////////////////////////////////////////////////////////
// File: midpoint_method.c                                                    //
// Routines:                                                                  //
//    Midpoint_Method                                                         //
////////////////////////////////////////////////////////////////////////////////

////////////////////////////////////////////////////////////////////////////////
//                                                                            //
//  Description:                                                              //
//     The midpoint method for approximating the solution of the differential //
//     equation y'(x) = f(x,y) with initial condition y(a) = c numerically    //
//     approximates the derivative y'(x) by ( y(x+h) - y(x-h) ) / 2h.         //
//                                                                            //
//     Let y[n] = y(a + nh), x[n] = a + nh, and f[n] = f(x[n],y[n]), then the //
//     recursion formula for the midpoint method is:                          //
//                          y[n+1] = y[n-1] + 2h * f[n].                      //
//                                                                            //
//     The midpoint method is a two step method and so requires two initial   //
//     values y[0] at x = a and y[1] at x = a + h to start the recursion.     //
//                                                                            //
//     Expand y(x+h) in a Taylor series about x,                              //
//       y(x+h) = y(x) + h*y'(x) + (h^2 / 2)*y''(x) + (h^3 / 6)*y'''(x) +...  //
//     and subtract y(x-h) from y(x+h)                                        //
//       y(x+h)-y(x-h) = 2h*y'(x) + 2(h^3 / 6)*y'''(x) + ... .                //
//                                                                            //
//     Thus the midpoint method is locally a third order and globally a       //
//     a second order procedure.                                              //
//                                                                            //
//     The midpoint method is preferrable to Euler's method, but as with      //
//     Euler's method it should not be used or even considered except in      //
//     cases for which only a few steps are desired.                          //
//                                                                            //
////////////////////////////////////////////////////////////////////////////////

////////////////////////////////////////////////////////////////////////////////
//  void Midpoint_Method( double (*f)(double, double), double y[], double x0, //
//                                          double h, int number_of_steps );  //
//                                                                            //
//  Description:                                                              //
//     This routine uses the midpoint method to approximate the solution of   //
//     the differential equation y'=f(x,y) with the initial condition y = y[0]//
//     at x = a and recursion starting value y = y[1] at x = a+h.  The values //
//     are returned in y[n], the value of y evaluated at x = a + n * h.       //
//                                                                            //
//  Arguments:                                                                //
//     double *f  Pointer to the integrand, a function of a two variables     //
//                of type double.                                             //
//     double y[] On input y[0] is the initial value of y at x = a, and y[1]  //
//                is the value of y at x = a + h, on output for i > 1         //
//                y[i] = y[i-2] + 2h f(x[i],y[i]), where x[i] = a + i * h.    //
//     double a   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) from x = a to x = a + number_of_steps * h  //
//     is stored in the input array y[].                                      //
//                                                                            //
////////////////////////////////////////////////////////////////////////////////
//                                                                            //
void Midpoint_Method( double (*f)(double, double), double y[], double a,
                                             double h, int number_of_steps ) {

   double x = a + h;
   double h2 = h + h;
   int i;

   for (i = 1; i < number_of_steps; x += h, i++) 
      y[i+1] = y[i-1] +  h2 * f( x, y[i] );
}