Quadrature: spec tweaks & Simpson's rule

doc/specs/stdlib_experimental_quadrature.md

@@ -76,26 +76,26 @@ program demo_trapz_weights
     real :: x(5) = [0., 1., 2., 3., 4.]
     real :: y(5) = x**2
     real :: w(5) 
-    w = trapz_weight(x)
-    print *, dot_product(w, y)
+    w = trapz_weights(x)
+    print *, sum(w*y)
 ! 22.0
 end program demo_trapz_weights
 
 ```
 
-## `simps` - integrate sampled values using Simpson's rule (to be implemented)
+## `simps` - integrate sampled values using Simpson's rule
 
 ### Description
 
 Returns the Simpson's rule integral of an array `y` representing discrete samples of a function. The integral is computed assuming either equidistant abscissas with spacing `dx` or arbitary abscissas `x`. 
 
-Simpson's rule is defined for odd-length arrays only. For even-length arrays, an optional argument `even` may be used to specify at which index to replace Simpson's rule with Simpson's 3/8 rule. The 3/8 rule will be used for the array section `y(even:even+4)` and the ordinary Simpon's rule will be used elsewhere.
+Simpson's ordinary ("1/3") rule is used for odd-length arrays. For even-length arrays, Simpson's 3/8 rule is also utilized in a way that depends on the value of `even`. If `even` is negative (positive), the 3/8 rule is used at the beginning (end) of the array. If `even` is zero or not present, the result is as if the 3/8 rule were first used at the beginning of the array, then at the end of the array, and these two results were averaged.
 
 ### Syntax
 
-`result = simps(y, x [, even])`
+`result = [[stdlib_experimental_quadrature(module):simps(interface)]](y, x [, even])`
 
-`result = simps(y, dx [, even])`
+`result = [[stdlib_experimental_quadrature(module):simps(interface)]](y, dx [, even])`
 
 ### Arguments
 
@@ -105,37 +105,48 @@ Simpson's rule is defined for odd-length arrays only. For even-length arrays, an
 
 `dx`: Shall be a scalar of type `real` having the same kind as `y`.
 
-`even`: (Optional) Shall be a scalar integer of default kind. Its default value is `1`.
+`even`: (Optional) Shall be a default-kind `integer`.
 
 ### Return value
 
 The result is a scalar of type `real` having the same kind as `y`.
 
 If the size of `y` is zero or one, the result is zero.
 
-If the size of `y` is two, the result is the same as if `trapz` had been called instead, regardless of the value of `even`.
+If the size of `y` is two, the result is the same as if `trapz` had been called instead.
 
 ### Example
 
-TBD
+```fortran
+program demo_simps
+    use stdlib_experimental_quadrature, only: simps
+    implicit none
+    real :: x(5) = [0., 1., 2., 3., 4.]
+    real :: y(5) = 3.*x**2
+    print *, simps(y, x) 
+! 64.0
+    print *, simps(y, 0.5) 
+! 32.0
+end program demo_simps
+```
 
-## `simps_weights` - Simpson's rule weights for given abscissas (to be implemented)
+## `simps_weights` - Simpson's rule weights for given abscissas
 
 ### Description
 
 Given an array of abscissas `x`, computes the array of weights `w` such that if `y` represented function values tabulated at `x`, then `sum(w*y)` produces a Simpson's rule approximation to the integral.
 
-Simpson's rule is defined for odd-length arrays only. For even-length arrays, an optional argument `even` may be used to specify at which index to replace Simpson's rule with Simpson's 3/8 rule. The 3/8 rule will be used for the array section `x(even:even+4)` and the ordinary Simpon's rule will be used elsewhere.
+Simpson's ordinary ("1/3") rule is used for odd-length arrays. For even-length arrays, Simpson's 3/8 rule is also utilized in a way that depends on the value of `even`. If `even` is negative (positive), the 3/8 rule is used at the beginning (end) of the array and the 1/3 rule used elsewhere. If `even` is zero or not present, the result is as if the 3/8 rule were first used at the beginning of the array, then at the end of the array, and then these two results were averaged.
 
 ### Syntax
 
-`result = simps_weights(x [, even])`
+`result = [[stdlib_experimental_quadrature(module):simps_weights(interface)]](x [, even])`
 
 ### Arguments
 
 `x`: Shall be a rank-one array of type `real`.
 
-`even`: (Optional) Shall be a scalar integer of default kind. Its default value is `1`.
+`even`: (Optional) Shall be a default-kind `integer`.
 
 ### Return value
 
@@ -147,4 +158,15 @@ If the size of `x` is two, then the result is the same as if `trapz_weights` had
 
 ### Example
 
-TBD
+```fortran
+program demo_simps_weights
+    use stdlib_experimental_quadrature, only: simps_weights
+    implicit none
+    real :: x(5) = [0., 1., 2., 3., 4.]
+    real :: y(5) = 3.*x**2
+    real :: w(5) 
+    w = simps_weights(x)
+    print *, sum(w*y)
+! 64.0
+end program demo_simps_weights
+```

src/CMakeLists.txt

@@ -13,6 +13,7 @@ set(fppFiles
     stdlib_experimental_stats_var.fypp
     stdlib_experimental_quadrature.fypp
     stdlib_experimental_quadrature_trapz.fypp
+    stdlib_experimental_quadrature_simps.fypp
 )
 
 

src/stdlib_experimental_quadrature.fypp

@@ -1,5 +1,4 @@
 #:include "common.fypp"
-
 module stdlib_experimental_quadrature
     !! ([Specification](../page/specs/stdlib_experimental_quadrature.html#description))
     use stdlib_experimental_kinds, only: sp, dp, qp
@@ -18,63 +17,82 @@ module stdlib_experimental_quadrature
     interface trapz
         !! Integrates sampled values using trapezoidal rule
         !! ([Specification](../page/specs/stdlib_experimental_quadrature.html#description))
-        #:for KIND in REAL_KINDS
-        pure module function trapz_dx_${KIND}$(y, dx) result(integral)
-            real(${KIND}$), dimension(:), intent(in) :: y
-            real(${KIND}$), intent(in) :: dx
-            real(${KIND}$) :: integral
-        end function trapz_dx_${KIND}$
+        #:for k1, t1 in REAL_KINDS_TYPES
+        pure module function trapz_dx_${k1}$(y, dx) result(integral)
+          ${t1}$, dimension(:), intent(in) :: y
+          ${t1}$, intent(in) :: dx
+          ${t1}$ :: integral
+        end function trapz_dx_${k1}$      
         #:endfor
-        #:for KIND in REAL_KINDS
-        pure module function trapz_x_${KIND}$(y, x) result(integral)
-            real(${KIND}$), dimension(:), intent(in) :: y
-            real(${KIND}$), dimension(size(y)), intent(in) :: x
-            real(${KIND}$) :: integral
-        end function trapz_x_${KIND}$
+        #:for k1, t1 in REAL_KINDS_TYPES
+        module function trapz_x_${k1}$(y, x) result(integral)
+            ${t1}$, dimension(:), intent(in) :: y
+            ${t1}$, dimension(:), intent(in) :: x
+            ${t1}$ :: integral
+        end function trapz_x_${k1}$
         #:endfor
     end interface trapz
 
 
     interface trapz_weights
         !! Integrates sampled values using trapezoidal rule weights for given abscissas
         !! ([Specification](../page/specs/stdlib_experimental_quadrature.html#description_1))
-        #:for KIND in REAL_KINDS
-        pure module function trapz_weights_${KIND}$(x) result(w)
-            real(${KIND}$), dimension(:), intent(in) :: x
-            real(${KIND}$), dimension(size(x)) :: w
-        end function trapz_weights_${KIND}$
+        #:for k1, t1 in REAL_KINDS_TYPES
+        pure module function trapz_weights_${k1}$(x) result(w)
+            ${t1}$, dimension(:), intent(in) :: x
+            ${t1}$, dimension(size(x)) :: w
+        end function trapz_weights_${k1}$
         #:endfor
     end interface trapz_weights
 
 
     interface simps
-        #:for KIND in REAL_KINDS
-        pure module function simps_dx_${KIND}$(y, dx, even) result(integral)
-            real(${KIND}$), dimension(:), intent(in) :: y
-            real(${KIND}$), intent(in) :: dx
+        !! Integrates sampled values using Simpson's rule
+        !! ([Specification](../page/specs/stdlib_experimental_quadrature.html#description_3))
+        ! "recursive" is an implementation detail
+        #: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
-            real(${KIND}$) :: integral
-        end function simps_dx_${KIND}$
+            ${t1}$ :: integral
+        end function simps_dx_${k1}$
         #:endfor
-        #:for KIND in REAL_KINDS
-        pure module function simps_x_${KIND}$(y, x, even) result(integral)
-            real(${KIND}$), dimension(:), intent(in) :: y
-            real(${KIND}$), dimension(size(y)), intent(in) :: x
+        #: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
-            real(${KIND}$) :: integral
-        end function simps_x_${KIND}$
+            ${t1}$ :: integral
+        end function simps_x_${k1}$
         #:endfor
     end interface simps
 
 
     interface simps_weights
-        #:for KIND in REAL_KINDS
-        pure module function simps_weights_${KIND}$(x, even) result(w)
-            real(${KIND}$), dimension(:), intent(in) :: x
-            real(${KIND}$), dimension(size(x)) :: w
+        !! Integrates sampled values using trapezoidal rule weights for given abscissas
+        !! ([Specification](../page/specs/stdlib_experimental_quadrature.html#description_3))
+        #: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
-        end function simps_weights_${KIND}$
+            ${t1}$, dimension(size(x)) :: w
+        end function simps_weights_${k1}$
         #:endfor
     end interface simps_weights
 
+
+    ! Interface for a simple f(x)-style integrand function.
+    ! Could become fancier as we learn about the performance
+    ! ramifications of different ways to do callbacks.
+    abstract interface
+        #:for k1, t1 in REAL_KINDS_TYPES
+        pure function integrand_${k1}$(x) result(f)
+            import :: ${k1}$
+            ${t1}$, intent(in) :: x
+            ${t1}$ :: f
+        end function integrand_${k1}$
+        #:endfor
+    end interface
+
 end module stdlib_experimental_quadrature

src/stdlib_experimental_quadrature_simps.fypp

@@ -0,0 +1,278 @@
+#:include "common.fypp"
+
+submodule (stdlib_experimental_quadrature) stdlib_experimental_quadrature_simps
+    use stdlib_experimental_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_experimental_quadrature_simps

src/stdlib_experimental_quadrature_trapz.fypp

@@ -1,7 +1,7 @@
 #:include "common.fypp"
 
 submodule (stdlib_experimental_quadrature) stdlib_experimental_quadrature_trapz
-
+    use stdlib_experimental_error, only: check
     implicit none
 
 contains
@@ -29,16 +29,17 @@ contains
 
     #:endfor
     #:for KIND in REAL_KINDS
-    
-    pure module function trapz_x_${KIND}$(y, x) result(integral)
+
+    module function trapz_x_${KIND}$(y, x) result(integral)
         real(${KIND}$), dimension(:), intent(in) :: y
-        real(${KIND}$), dimension(size(y)), intent(in) :: x
+        real(${KIND}$), dimension(:), intent(in) :: x
         real(${KIND}$) :: integral
 
         integer :: i
         integer :: n
 
         n = size(y)
+        call check(size(x) == n, "trapz: Arguments `x` and `y` must be the same size.")
 
         select case (n)
         case (0:1)

src/tests/quadrature/CMakeLists.txt

@@ -1,2 +1,2 @@
 ADDTEST(trapz)
-
+ADDTEST(simps)