stdlib_quadrature_simps.fypp Source File


This file depends on

sourcefile~~stdlib_quadrature_simps.fypp~~EfferentGraph sourcefile~stdlib_quadrature_simps.fypp stdlib_quadrature_simps.fypp sourcefile~stdlib_quadrature.fypp stdlib_quadrature.fypp sourcefile~stdlib_quadrature_simps.fypp->sourcefile~stdlib_quadrature.fypp sourcefile~stdlib_error.f90 stdlib_error.f90 sourcefile~stdlib_quadrature_simps.fypp->sourcefile~stdlib_error.f90 sourcefile~stdlib_kinds.f90 stdlib_kinds.f90 sourcefile~stdlib_quadrature.fypp->sourcefile~stdlib_kinds.f90 sourcefile~stdlib_optval.fypp stdlib_optval.fypp sourcefile~stdlib_error.f90->sourcefile~stdlib_optval.fypp sourcefile~stdlib_optval.fypp->sourcefile~stdlib_kinds.f90

Contents


Source Code

#:include "common.fypp"

submodule (stdlib_quadrature) stdlib_quadrature_simps
    use stdlib_error, only: check
    implicit none

    ! internal use only
    interface simps38
    #:for k1, t1 in REAL_KINDS_TYPES
        module procedure simps38_dx_${k1}$
        module procedure simps38_x_${k1}$
    #:endfor
    end interface simps38

    ! internal use only
    interface simps38_weights
    #:for k1, t1 in REAL_KINDS_TYPES
        module procedure simps38_weights_${k1}$
    #:endfor
    end interface simps38_weights

contains

#:for k1, t1 in REAL_KINDS_TYPES

    pure recursive module function simps_dx_${k1}$(y, dx, even) result(integral)
        ${t1}$, dimension(:), intent(in) :: y
        ${t1}$, intent(in) :: dx
        integer, intent(in), optional :: even
        ${t1}$ :: integral

        integer :: n

        n = size(y)

        select case (n)
        case (0:1)
            integral = 0.0_${k1}$
        case (2)
            integral = 0.5_${k1}$*dx*(y(1) + y(2))
        case (3)
            integral = dx/3.0_${k1}$*(y(1) + 4*y(2) + y(3))
        case (4)
            integral = simps38(y, dx)
        ! case (5) not needed; handled by default
        case (6) ! needs special handling because of averaged 3/8's rule case
            if (present(even)) then
                if (even < 0) then
                    ! 3/8 rule on left
                    integral =  simps38(y(1:4), dx) + simps(y(4:6), dx)
                    return
                else if (even > 0) then
                    ! 3/8 rule on right
                    integral = simps(y(1:3), dx) + simps38(y(3:6), dx)
                    return
                else
                    ! fall through
                end if
            end if
            ! either `even` not present or is zero
            ! equivalent to averaging left and right
            integral = dx/48.0_${k1}$ * (17*(y(1) + y(6)) + 59*(y(2) + y(5)) + 44*(y(3) + y(4)))
        case default
            if (mod(n, 2) == 1) then
                integral = dx/3.0_${k1}$*(y(1) + 4*sum(y(2:n-1:2)) + 2*sum(y(3:n-2:2)) + y(n))
            else
                if (present(even)) then
                    if (even < 0) then
                        ! 3/8th rule on left
                        integral = simps38(y(1:4), dx) + simps(y(4:n), dx)
                        return
                    else if (even > 0) then
                        ! 3/8 rule on right
                        integral = simps(y(1:n-3), dx) + simps38(y(n-3:n), dx)
                        return
                    else
                        ! fall through
                    end if
                end if
                ! either `even` not present or is zero
                ! equivalent to averaging left and right
                integral = dx/48.0_${k1}$ * (17*(y(1) + y(n)) + 59*(y(2) + y(n-1)) &
                        + 43*(y(3) + y(n-2)) + 49*(y(4) + y(n-3)) + 48*sum(y(5:n-4)))
            end if
        end select
    end function simps_dx_${k1}$

#:endfor
#:for k1, t1 in REAL_KINDS_TYPES

    recursive module function simps_x_${k1}$(y, x, even) result(integral)
        ${t1}$, dimension(:), intent(in) :: y
        ${t1}$, dimension(:), intent(in) :: x
        integer, intent(in), optional :: even
        ${t1}$ :: integral

        integer :: i
        integer :: n

        ${t1}$ :: h1, h2
        ${t1}$ :: a, b, c

        n = size(y)
        call check(size(x) == n, "simps: Arguments `x` and `y` must be the same size.")

        select case (n)
        case (0:1)
            integral = 0.0_${k1}$
        case (2)
            integral = 0.5_${k1}$*(x(2) - x(1))*(y(1) + y(2))
        case (3)
            h1 = x(2) - x(1)
            h2 = x(3) - x(2)
            a = (2*h1**2 + h1*h2 - h2**2)/(6*h1)
            b = (h1+h2)**3/(6*h1*h2)
            c = (2*h2**2 + h1*h2 - h1**2)/(6*h2)
            integral = a*y(1) + b*y(2) + c*y(3)
        case (4)
            integral = simps38(y, x)
        ! case (6) unneeded; handled by default
        case default
            if (mod(n, 2) == 1) then
                integral = 0.0_${k1}$
                do i = 1, n-2, 2
                    h1 = x(i+1) - x(i)
                    h2 = x(i+2) - x(i+1)
                    a = (2*h1**2 + h1*h2 - h2**2)/(6*h1)
                    b = (h1+h2)**3/(6*h1*h2)
                    c = (2*h2**2 + h1*h2 - h1**2)/(6*h2)
                    integral = integral + a*y(i) + b*y(i+1) + c*y(i+2)
                end do
            else
                if (present(even)) then
                    if (even < 0) then
                        ! 3/8 rule on left
                        integral = simps38(y(1:4), x(1:4)) + simps(y(4:n), x(4:n))
                        return
                    else if (even > 0) then
                        ! 3/8 rule on right
                        integral = simps(y(1:n-3), x(1:n-3)) + simps38(y(n-3:n), x(n-3:n))
                        return
                    else
                        ! fall through
                    end if
                end if
                ! either `even` not present or is zero
                integral = 0.5_${k1}$ * ( simps38(y(1:4), x(1:4)) + simps(y(4:n), x(4:n)) &
                    + simps(y(1:n-3), x(1:n-3)) + simps38(y(n-3:n), x(n-3:n)) )
            end if
        end select
    end function simps_x_${k1}$

#:endfor
#:for k1, t1 in REAL_KINDS_TYPES

    pure recursive module function simps_weights_${k1}$(x, even) result(w)
        ${t1}$, dimension(:), intent(in) :: x
        integer, intent(in), optional :: even
        ${t1}$, dimension(size(x)) :: w

        integer :: i, n
        ${t1}$ :: h1, h2

        n = size(x)

        select case (n)
        case (0)
            ! no action needed
        case (1)
            w(1) = 0.0_${k1}$
        case (2)
            w = 0.5_${k1}$*(x(2) - x(1))
        case (3)
            h1 = x(2) - x(1)
            h2 = x(3) - x(2)
            w(1) = (2*h1**2 + h1*h2 - h2**2)/(6*h1)
            w(2) = (h1+h2)**3/(6*h1*h2)
            w(3) = (2*h2**2 + h1*h2 - h1**2)/(6*h2)
        case (4)
            w = simps38_weights(x)
        case default
            if (mod(n, 2) == 1) then
                w = 0.0_${k1}$
                do i = 1, n-2, 2
                    h1 = x(i+1) - x(i)
                    h2 = x(i+2) - x(i+1)
                    w(i) = w(i) + (2*h1**2 + h1*h2 - h2**2)/(6*h1)
                    w(i+1) = w(i+1) + (h1+h2)**3/(6*h1*h2)
                    w(i+2) = w(i+2) + (2*h2**2 + h1*h2 - h1**2)/(6*h2)
                end do
            else
                if (present(even)) then
                    if (even < 0) then
                        ! 3/8 rule on left
                        w = 0.0_${k1}$
                        w(1:4) = simps38_weights(x(1:4))
                        w(4:n) = w(4:n) + simps_weights(x(4:n)) ! position 4 needs both rules
                        return
                    else if (even > 0) then
                        ! 3/8 rule on right
                        w = 0.0_${k1}$
                        w(1:n-3) = simps_weights(x(1:n-3))
                        w(n-3:n) = w(n-3:n) + simps38_weights(x(n-3:n)) ! position n-3 needs both rules
                        return
                    else
                        ! fall through
                    end if
                end if
                ! either `even` not present or is zero
                w = 0.0_${k1}$
                ! 3/8 rule on left
                w(1:4) = simps38_weights(x(1:4))
                w(4:n) = w(4:n) + simps_weights(x(4:n))
                ! 3/8 rule on right
                w(1:n-3) = w(1:n-3) + simps_weights(x(1:n-3))
                w(n-3:n) = w(n-3:n) + simps38_weights(x(n-3:n))
                ! average
                w = 0.5_${k1}$ * w
            end if
        end select
    end function simps_weights_${k1}$

#:endfor
#:for k1, t1 in REAL_KINDS_TYPES

    pure function simps38_dx_${k1}$(y, dx) result (integral)
        ${t1}$, dimension(4), intent(in) :: y
        ${t1}$, intent(in) :: dx
        ${t1}$ :: integral

        integral = 3.0_${k1}$*dx/8.0_${k1}$ * (y(1) + y(4) + 3*(y(2) + y(3)))
    end function simps38_dx_${k1}$

#:endfor
#: for k1, t1 in REAL_KINDS_TYPES

    pure function simps38_x_${k1}$(y, x) result(integral)
        ${t1}$, dimension(4), intent(in) :: y
        ${t1}$, dimension(4), intent(in) :: x
        ${t1}$ :: integral

        ${t1}$ :: h1, h2, h3
        ${t1}$ :: a, b, c, d

        h1 = x(2) - x(1)
        h2 = x(3) - x(2)
        h3 = x(4) - x(3)

        a = (h1+h2+h3)*(3*h1**2 + 2*h1*h2 - 2*h1*h3 - h2**2 + h3**2)/(12*h1*(h1+h2))
        b = (h1+h2-h3)*(h1+h2+h3)**3/(12*h1*h2*(h2+h3))
        c = (h2+h3-h1)*(h1+h2+h3)**3/(12*h2*h3*(h1+h2))
        d = (h1+h2+h3)*(3*h3**2 + 2*h2*h3 - 2*h1*h3 - h2**2 + h1**2)/(12*h3*(h2+h3))

        integral = a*y(1) + b*y(2) + c*y(3) + d*y(4)
    end function simps38_x_${k1}$

#:endfor
#:for k1, t1 in REAL_KINDS_TYPES

    pure function simps38_weights_${k1}$(x) result(w)
        ${t1}$, intent(in) :: x(4)
        ${t1}$ :: w(size(x))

        ${t1}$ :: h1, h2, h3

        h1 = x(2) - x(1)
        h2 = x(3) - x(2)
        h3 = x(4) - x(3)

        w(1) = (h1+h2+h3)*(3*h1**2 + 2*h1*h2 - 2*h1*h3 - h2**2 + h3**2)/(12*h1*(h1+h2))
        w(2) = (h1+h2-h3)*(h1+h2+h3)**3/(12*h1*h2*(h2+h3))
        w(3) = (h2+h3-h1)*(h1+h2+h3)**3/(12*h2*h3*(h1+h2))
        w(4) = (h1+h2+h3)*(3*h3**2 + 2*h2*h3 - 2*h1*h3 - h2**2 + h1**2)/(12*h3*(h2+h3))
    end function simps38_weights_${k1}$

#:endfor

end submodule stdlib_quadrature_simps