Merge pull request #146 from nshaffer/dev-quadrature

First steps toward quadrature

Author: milancurcic

src/CMakeLists.txt

@@ -7,6 +7,8 @@ set(fppFiles
     stdlib_experimental_stats.fypp
     stdlib_experimental_stats_mean.fypp
     stdlib_experimental_stats_var.fypp
+    stdlib_experimental_quadrature.fypp
+    stdlib_experimental_quadrature_trapz.fypp
 )
 
 

src/stdlib_experimental_quadrature.fypp

@@ -0,0 +1,75 @@
+#:include "common.fypp"
+
+module stdlib_experimental_quadrature
+    use stdlib_experimental_kinds, only: sp, dp, qp
+
+    implicit none
+
+    private
+
+    ! array integration
+    public :: trapz
+    public :: trapz_weights
+    public :: simps
+    public :: simps_weights
+
+
+    interface trapz
+        #: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}$
+        #: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}$
+        #:endfor
+    end interface trapz
+
+
+    interface trapz_weights
+        #: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}$
+        #: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
+            integer, intent(in), optional :: even
+            real(${KIND}$) :: integral
+        end function simps_dx_${KIND}$
+        #: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
+            integer, intent(in), optional :: even
+            real(${KIND}$) :: integral
+        end function simps_x_${KIND}$
+        #: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
+            integer, intent(in), optional :: even
+        end function simps_weights_${KIND}$
+        #:endfor
+    end interface simps_weights
+
+end module stdlib_experimental_quadrature

src/stdlib_experimental_quadrature.md

@@ -0,0 +1,141 @@
+# Numerical integration
+
+## Implemented
+
+* `trapz`
+* `trapz_weights`
+
+## `trapz` - integrate sampled values using trapezoidal rule
+
+Returns the trapezoidal 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`.
+
+### Syntax
+
+`result = trapz(y, x)`
+
+`result = trapz(y, dx)`
+
+### Arguments
+
+`y`: Shall be a rank-one array of type `real`.
+
+`x`: Shall be a rank-one array of type `real` having the same kind and size as `y`. 
+
+`dx`: Shall be a scalar of type `real` having the same kind as `y`.
+
+### 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.
+
+### Example
+
+```fortran
+program demo_trapz
+    use stdlib_experimental_quadrature, only: trapz
+    implicit none
+    real :: x(5) = [0., 1., 2., 3., 4.]
+    real :: y(5) = x**2
+    print *, trapz(y, x) 
+! 22.0
+    print *, trapz(y, 0.5) 
+! 11.0
+end program
+```
+
+## `trapz_weights` - trapezoidal rule weights for given abscissas
+
+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 trapezoidal rule approximation to the integral.
+
+### Syntax
+
+`result = trapz_weights(x)`
+
+### Arguments
+
+`x`: Shall be a rank-one array of type `real`.
+
+### Return value
+
+The result is a `real` array with the same size and kind as `x`.
+
+If the size of `x` is one, then the sole element of the result is zero.
+
+### Example
+
+```fortran
+program demo_trapz_weights
+    use stdlib_experimental_quadrature, only: trapz_weights
+    implicit none
+    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)
+! 22.0
+end program
+
+```
+
+# `simps` - integrate sampled values using Simpson's rule
+
+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.
+
+### Syntax
+
+`result = simps(y, x [, even])`
+
+`result = simps(y, dx [, even])`
+
+### Arguments
+
+`y`: Shall be a rank-one array of type `real`.
+
+`x`: Shall be a rank-one array of type `real` having the same kind and size as `y`. 
+
+`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`.
+
+### 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`.
+
+### Example
+
+TBD
+
+# `simps_weights` - Simpson's rule weights for given abscissas
+
+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.
+
+### Syntax
+
+`result = simps_weights(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`.
+
+### Return value
+
+The result is a `real` array with the same size and kind as `x`.
+
+If the size of `x` is one, then the sole element of the result is zero.
+
+If the size of `x` is two, then the result is the same as if `trapz_weights` had been called instead.
+
+### Example
+
+TBD

src/stdlib_experimental_quadrature_trapz.fypp

@@ -0,0 +1,86 @@
+#:include "common.fypp"
+
+submodule (stdlib_experimental_quadrature) stdlib_experimental_quadrature_trapz
+
+    implicit none
+
+contains
+
+    #: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
+
+        integer :: n
+
+        n = size(y)
+        
+        select case (n)
+        case (0:1)
+            integral = 0.0_${KIND}$
+        case (2)
+            integral = 0.5_${KIND}$*dx*(y(1) + y(2))
+        case default
+            integral = dx*(sum(y(2:n-1)) + 0.5_${KIND}$*(y(1) + y(n)))
+        end select
+    end function trapz_dx_${KIND}$
+    
+    #: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
+
+        integer :: i
+        integer :: n
+
+        n = size(y)
+
+        select case (n)
+        case (0:1)
+            integral = 0.0_${KIND}$
+        case (2)
+            integral = 0.5_${KIND}$*(x(2) - x(1))*(y(1) + y(2))
+        case default
+            integral = 0.0_${KIND}$
+            do i = 2, n
+                integral = integral + (x(i) - x(i-1))*(y(i) + y(i-1))
+            end do
+            integral = 0.5_${KIND}$*integral
+        end select
+    end function trapz_x_${KIND}$
+
+    #:endfor
+    #: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
+
+        integer :: i
+        integer :: n
+
+        n = size(x)
+
+        select case (n)
+        case (0)
+            ! no action needed
+        case (1)
+            w(1) = 0.0_${KIND}$
+        case (2)
+            w = 0.5_${KIND}$*(x(2) - x(1))
+        case default
+            w(1) = 0.5_${KIND}$*(x(2) - x(1))
+            w(n) = 0.5_${KIND}$*(x(n) - x(n-1))
+            do i = 2, size(x)-1
+                w(i) = 0.5_${KIND}$*(x(i+1) - x(i-1))
+            end do
+        end select
+    end function trapz_weights_${KIND}$
+
+#:endfor
+end submodule stdlib_experimental_quadrature_trapz

src/tests/CMakeLists.txt

@@ -11,6 +11,7 @@ add_subdirectory(io)
 add_subdirectory(optval)
 add_subdirectory(stats)
 add_subdirectory(system)
+add_subdirectory(quadrature)
 
 ADDTEST(always_skip)
 set_tests_properties(always_skip PROPERTIES SKIP_RETURN_CODE 77)

src/tests/quadrature/CMakeLists.txt

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

src/tests/quadrature/Makefile.manual

@@ -0,0 +1,3 @@
+PROGS_SRC = test_trapz.f90
+
+include ../Makefile.manual.test.mk