The Modified Midpoint Method

The Modified Midpoint Method is based on the following formulas:

$$\begin{eqnarray*}x(t_0) &=& x_0 \\ x_1 = x(t_1) &=& x_0 + f(x_0, t_0) \Delta t... ..._{m-1} + f(x_m, t_m) 2\Delta t \quad t_m = t_{m-1} + \Delta t \end{eqnarray*}$$

This method is also called a leap-frog method and is often used in particle codes, being quite a lot cheaper than the Runge-Kutta method, while at the same time, offering a better stability than the Euler method.

The following code illustrates Fortran-90 implementation:

SUBROUTINE mmid(y, dydx, xs, htot, nstep, yout, derivs)

  USE nrtype; USE nrutil, ONLY: assert_eq, swap
  IMPLICIT NONE
  INTEGER(i4b), INTENT(in) :: nstep
  REAL(sp), INTENT(in) :: xs, htot
  REAL(sp), DIMENSION(:), INTENT(in) :: y, dydx
  REAL(sp), DIMENSION(:), INTENT(out) :: yout
  INTERFACE
     SUBROUTINE derivs(x, y, dydx)
       USE nrtype
       IMPLICIT NONE
       REAL(sp), INTENT(in) :: x
       REAL(sp), DIMENSION(:), INTENT(in) :: y
       REAL(sp), DIMENSION(:), INTENT(out) :: dydx
     END SUBROUTINE derivs
  END INTERFACE
  INTEGER(i4b) :: n, ndum
  REAL(sp) :: h, h2, x
  REAL(sp), DIMENSION(SIZE(y)) :: ym, yn

  ndum = assert_eq(SIZE(y), SIZE(dydx), SIZE(yout), 'mmid')
  h = htot / nstep
  ym = y
  yn = y + h * dydx
  x = xs + h
  CALL derivs(x, yn, yout)
  h2 = 2.0_sp * h
  DO n = 2, nstep
     CALL swap(ym, yn)
     yn = yn + h2 * yout
     x = x + h
     CALL derivs(x, yn, yout)
  END DO
  yout = 0.5_sp * (ym + yn + h * yout)

END SUBROUTINE mmid