Fortran Wiki
integration

PCHIP method (accurate)

Below is a function that (i) interpolates an array using a Piecewise Cubic Hermitian Interpolator (using DPCHEZ), and then (ii) integrates the interpolation (using DPCHQA). The required external methods are provided by the PCHIP module.

This method is quite accurate, but not as fast as e.g. the trapezoid method below.

  function integrate(x, y, a, b) result(r)
    !! This function constructs a piecewise cubic Hermitian interpolation of an array y(x) based on
    !! discrete numerical data, and subsequently evaluates the integral of the interpolation in the
    !! range (a,b). Note that the mesh spacing of x does not necessarily have to be uniform.
    real(wp), intent(in)  :: x(:)        !! Variable x
    real(wp), intent(in)  :: y(size(x))  !! Function y(x)
    real(wp), intent(in)  :: a           !! Left endpoint
    real(wp), intent(in)  :: b           !! Right endpoint
    real(wp)              :: r           !! Integral ∫y(x)·dx

    real(wp), external    :: dpchqa
    real(wp)              :: d(size(x))
    integer               :: err

    ! Create a PCHIP interpolation of the input data
    call dpchez(size(x), x, y, d, .false., 0, 0, err)

    ! Integrate the interpolation in the provided range
    r = dpchqa(size(x), x, y, d, a, b, err)
  end function

Trapezoid method (fast)

Below is a simple function for numerically calculating the integral of an an array using the trapezoid method. Since it is expressed using whole-array operations, a good compiler should be able to vectorize it automatically, rendering it very fast.

  pure function integrate(x, y) result(r)
    !! Calculates the integral of an array y with respect to x using the trapezoid
    !! approximation. Note that the mesh spacing of x does not have to be uniform.
    real(wp), intent(in)  :: x(:)         !! Variable x
    real(wp), intent(in)  :: y(size(x))   !! Function y(x)
    real(wp)              :: r            !! Integral ∫y(x)·dx

    ! Integrate using the trapezoidal rule
    associate(n => size(x))
      r = sum((y(1+1:n-0) + y(1+0:n-1))*(x(1+1:n-0) - x(1+0:n-1)))/2
    end associate
  end function