From 97ce9d30aaf5f0b174df9ffb0f83f4453bc3df6f Mon Sep 17 00:00:00 2001
From: Nathaniel Shaffer <nrshaffer@protonmail.com>
Date: Thu, 30 Jan 2020 00:11:37 -0700
Subject: [PATCH 01/10] initial commit -- trapz
---
src/CMakeLists.txt | 2 +
src/stdlib_experimental_quadrature.fypp | 37 +++
src/stdlib_experimental_quadrature_trapz.fypp | 87 +++++++
src/tests/CMakeLists.txt | 1 +
src/tests/quadrature/CMakeLists.txt | 2 +
src/tests/quadrature/Makefile.manual | 3 +
src/tests/quadrature/test_trapz.f90 | 236 ++++++++++++++++++
7 files changed, 368 insertions(+)
create mode 100644 src/stdlib_experimental_quadrature.fypp
create mode 100644 src/stdlib_experimental_quadrature_trapz.fypp
create mode 100644 src/tests/quadrature/CMakeLists.txt
create mode 100644 src/tests/quadrature/Makefile.manual
create mode 100644 src/tests/quadrature/test_trapz.f90
diff --git a/src/CMakeLists.txt b/src/CMakeLists.txt
index 90bd8e1ed..480e2737a 100644
--- a/src/CMakeLists.txt
+++ b/src/CMakeLists.txt
@@ -5,6 +5,8 @@ set(fppFiles
stdlib_experimental_io.fypp
stdlib_experimental_stats.fypp
stdlib_experimental_stats_mean.fypp
+ stdlib_experimental_quadrature.fypp
+ stdlib_experimental_quadrature_trapz.fypp
)
diff --git a/src/stdlib_experimental_quadrature.fypp b/src/stdlib_experimental_quadrature.fypp
new file mode 100644
index 000000000..e049eaab9
--- /dev/null
+++ b/src/stdlib_experimental_quadrature.fypp
@@ -0,0 +1,37 @@
+#:include "common.fypp"
+
+module stdlib_experimental_quadrature
+ use stdlib_experimental_kinds, only: sp, dp, qp
+
+ implicit none
+
+ private
+ public :: trapz, trapz_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
+
+end module stdlib_experimental_quadrature
diff --git a/src/stdlib_experimental_quadrature_trapz.fypp b/src/stdlib_experimental_quadrature_trapz.fypp
new file mode 100644
index 000000000..8bdb91cbe
--- /dev/null
+++ b/src/stdlib_experimental_quadrature_trapz.fypp
@@ -0,0 +1,87 @@
+#:include "common.fypp"
+
+submodule (stdlib_experimental_quadrature) stdlib_experimental_quadrature_trapz
+ use stdlib_experimental_kinds, only: sp, dp, qp
+
+ 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
diff --git a/src/tests/CMakeLists.txt b/src/tests/CMakeLists.txt
index df5bd0a09..f3b7d434c 100644
--- a/src/tests/CMakeLists.txt
+++ b/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)
diff --git a/src/tests/quadrature/CMakeLists.txt b/src/tests/quadrature/CMakeLists.txt
new file mode 100644
index 000000000..08f1ce4c6
--- /dev/null
+++ b/src/tests/quadrature/CMakeLists.txt
@@ -0,0 +1,2 @@
+ADDTEST(trapz)
+
diff --git a/src/tests/quadrature/Makefile.manual b/src/tests/quadrature/Makefile.manual
new file mode 100644
index 000000000..83061a7bf
--- /dev/null
+++ b/src/tests/quadrature/Makefile.manual
@@ -0,0 +1,3 @@
+PROGS_SRC = test_trapz.f90
+
+include ../Makefile.manual.test.mk
diff --git a/src/tests/quadrature/test_trapz.f90 b/src/tests/quadrature/test_trapz.f90
new file mode 100644
index 000000000..b275643e8
--- /dev/null
+++ b/src/tests/quadrature/test_trapz.f90
@@ -0,0 +1,236 @@
+program test_trapz
+ use stdlib_experimental_kinds, only: sp, dp, qp
+ use stdlib_experimental_error, only: assert
+ use stdlib_experimental_quadrature, only: trapz, trapz_weights
+
+ implicit none
+
+ call test_trapz_sp
+ call test_trapz_dp
+ call test_trapz_qp
+
+ call test_trapz_weights_sp
+ call test_trapz_weights_dp
+ call test_trapz_weights_qp
+
+ call test_trapz_zero_sp
+ call test_trapz_zero_dp
+ call test_trapz_zero_qp
+
+contains
+
+ subroutine test_trapz_sp
+ integer, parameter :: n = 17
+ real(sp), dimension(n) :: y
+ real(sp), dimension(n) :: x
+ real(sp) :: val
+ real(sp) :: ans
+ integer :: i
+
+ print *, "test_trapz_sp"
+
+ y = [(real(i-1, sp), i = 1, n)]
+
+ val = trapz(y, 1.0_sp)
+ ans = 128.0_sp
+ call assert(val == ans)
+
+ val = trapz(y, 0.5_sp)
+ ans = 64.0_sp
+ call assert(val == ans)
+
+ x = [((i-1)*4.0_sp/real(n-1, sp), i = 1, n)]
+ val = trapz(y, x)
+ ans = 32.0_sp
+ call assert(val == ans)
+
+ x = y**2
+ val = trapz(y, x)
+ ans = 2728.0_sp
+ call assert(val == ans)
+ end subroutine test_trapz_sp
+
+ subroutine test_trapz_dp
+ integer, parameter :: n = 17
+ real(dp), dimension(n) :: y
+ real(dp), dimension(n) :: x
+ real(dp) :: val
+ real(dp) :: ans
+ integer :: i
+
+ print *, "test_trapz_dp"
+
+ y = [(real(i-1, dp), i = 1, n)]
+
+ val = trapz(y, 1.0_dp)
+ ans = 128.0_dp
+ call assert(val == ans)
+
+ val = trapz(y, 0.5_dp)
+ ans = 64.0_dp
+ call assert(val == ans)
+
+ x = [((i-1)*4.0_dp/real(n-1, dp), i = 1, n)]
+ val = trapz(y, x)
+ ans = 32.0_dp
+ call assert(val == ans)
+
+ x = y**2
+ val = trapz(y, x)
+ ans = 2728.0_sp
+ call assert(val == ans)
+ end subroutine test_trapz_dp
+
+
+ subroutine test_trapz_qp
+ integer, parameter :: n = 17
+ real(qp), dimension(n) :: y
+ real(qp), dimension(n) :: x
+ real(qp) :: val
+ real(qp) :: ans
+ integer :: i
+
+ print *, "test_trapz_qp"
+
+ y = [(real(i-1, qp), i = 1, n)]
+
+ val = trapz(y, 1.0_qp)
+ ans = 128.0_qp
+ call assert(val == ans)
+
+ val = trapz(y, 0.5_qp)
+ ans = 64.0_qp
+ call assert(val == ans)
+
+ x = [((i-1)*4.0_qp/real(n-1, qp), i = 1, n)]
+ val = trapz(y, x)
+ ans = 32.0_qp
+ call assert(val == ans)
+
+ x = y**2
+ val = trapz(y, x)
+ ans = 2728.0_qp
+ call assert(val == ans)
+ end subroutine test_trapz_qp
+
+
+ subroutine test_trapz_weights_sp
+ integer, parameter :: n = 17
+ real(sp), dimension(n) :: y
+ real(sp), dimension(n) :: x
+ real(sp), dimension(n) :: w
+ integer :: i
+ real(sp) :: val
+ real(sp) :: ans
+
+ print *, "test_trapz_weights_sp"
+
+ y = [(real(i-1, sp), i = 1, n)]
+
+ x = y
+ w = trapz_weights(x)
+ val = dot_product(w, y)
+ ans = trapz(y, x)
+ call assert(val == ans)
+
+ x = y**2
+ w = trapz_weights(x)
+ val = dot_product(w, y)
+ ans = trapz(y, x)
+ call assert(val == ans)
+
+ end subroutine test_trapz_weights_sp
+
+
+ subroutine test_trapz_weights_dp
+ integer, parameter :: n = 17
+ real(dp), dimension(n) :: y
+ real(dp), dimension(n) :: x
+ real(dp), dimension(n) :: w
+ integer :: i
+ real(dp) :: val
+ real(dp) :: ans
+
+ print *, "test_trapz_weights_dp"
+
+ y = [(real(i-1, dp), i = 1, n)]
+
+ x = y
+ w = trapz_weights(x)
+ val = dot_product(w, y)
+ ans = trapz(y, x)
+ call assert(val == ans)
+
+ x = y**2
+ w = trapz_weights(x)
+ val = dot_product(w, y)
+ ans = trapz(y, x)
+ call assert(val == ans)
+
+ end subroutine test_trapz_weights_dp
+
+
+ subroutine test_trapz_weights_qp
+ integer, parameter :: n = 17
+ real(qp), dimension(n) :: y
+ real(qp), dimension(n) :: x
+ real(qp), dimension(n) :: w
+ integer :: i
+ real(qp) :: val
+ real(qp) :: ans
+
+ print *, "test_trapz_weights_qp"
+
+ y = [(real(i-1, qp), i = 1, n)]
+
+ x = y
+ w = trapz_weights(x)
+ val = dot_product(w, y)
+ ans = trapz(y, x)
+ call assert(val == ans)
+
+ x = y**2
+ w = trapz_weights(x)
+ val = dot_product(w, y)
+ ans = trapz(y, x)
+ call assert(val == ans)
+
+ end subroutine test_trapz_weights_qp
+
+
+ subroutine test_trapz_zero_sp
+ real(sp), dimension(0) :: a
+
+ print *, "test_trapz_zero_sp"
+
+ call assert(trapz(a, 1.0_sp) == 0.0_sp)
+ call assert(trapz([1.0_sp], 1.0_sp) == 0.0_sp)
+ call assert(trapz(a, a) == 0.0_sp)
+ call assert(trapz([1.0_sp], [1.0_sp]) == 0.0_sp)
+ end subroutine test_trapz_zero_sp
+
+
+ subroutine test_trapz_zero_dp
+ real(dp), dimension(0) :: a
+
+ print *, "test_trapz_zero_dp"
+
+ call assert(trapz(a, 1.0_dp) == 0.0_dp)
+ call assert(trapz([1.0_dp], 1.0_dp) == 0.0_dp)
+ call assert(trapz(a, a) == 0.0_dp)
+ call assert(trapz([1.0_dp], [1.0_dp]) == 0.0_dp)
+ end subroutine test_trapz_zero_dp
+
+
+ subroutine test_trapz_zero_qp
+ real(qp), dimension(0) :: a
+
+ print *, "test_trapz_zero_qp"
+
+ call assert(trapz(a, 1.0_qp) == 0.0_qp)
+ call assert(trapz([1.0_qp], 1.0_qp) == 0.0_qp)
+ call assert(trapz(a, a) == 0.0_qp)
+ call assert(trapz([1.0_qp], [1.0_qp]) == 0.0_qp)
+ end subroutine test_trapz_zero_qp
+
+end program test_trapz
From 77a31aebb97e7dec5b5ca4a9e3cf088f37bfe037 Mon Sep 17 00:00:00 2001
From: Nathaniel Shaffer <nrshaffer@protonmail.com>
Date: Thu, 30 Jan 2020 00:18:28 -0700
Subject: [PATCH 02/10] bugfix? redundant use stmt
---
src/stdlib_experimental_quadrature_trapz.fypp | 1 -
1 file changed, 1 deletion(-)
diff --git a/src/stdlib_experimental_quadrature_trapz.fypp b/src/stdlib_experimental_quadrature_trapz.fypp
index 8bdb91cbe..0ea66fe5e 100644
--- a/src/stdlib_experimental_quadrature_trapz.fypp
+++ b/src/stdlib_experimental_quadrature_trapz.fypp
@@ -1,7 +1,6 @@
#:include "common.fypp"
submodule (stdlib_experimental_quadrature) stdlib_experimental_quadrature_trapz
- use stdlib_experimental_kinds, only: sp, dp, qp
implicit none
From 52386a6a95ecf3edc2dc9c9a6c1b62aaa64a6bcc Mon Sep 17 00:00:00 2001
From: Nathaniel Shaffer <nrshaffer@protonmail.com>
Date: Sun, 2 Feb 2020 08:26:56 -0700
Subject: [PATCH 03/10] Add simpson's rule API (not implemented)
---
src/stdlib_experimental_quadrature.fypp | 40 ++++++++++++++++++++++++-
1 file changed, 39 insertions(+), 1 deletion(-)
diff --git a/src/stdlib_experimental_quadrature.fypp b/src/stdlib_experimental_quadrature.fypp
index e049eaab9..ecb15840f 100644
--- a/src/stdlib_experimental_quadrature.fypp
+++ b/src/stdlib_experimental_quadrature.fypp
@@ -6,7 +6,13 @@ module stdlib_experimental_quadrature
implicit none
private
- public :: trapz, trapz_weights
+
+ ! array integration
+ public :: trapz
+ public :: trapz_weights
+ public :: simps
+ public :: simps_weights
+
interface trapz
#:for KIND in REAL_KINDS
@@ -25,6 +31,7 @@ module stdlib_experimental_quadrature
#:endfor
end interface trapz
+
interface trapz_weights
#:for KIND in REAL_KINDS
pure module function trapz_weights_${KIND}$(x) result(w)
@@ -34,4 +41,35 @@ module stdlib_experimental_quadrature
#: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
From b1fb2651ea1b04586baa0ee927ceb2c8fc04080a Mon Sep 17 00:00:00 2001
From: Nathaniel Shaffer <nrshaffer@protonmail.com>
Date: Wed, 5 Feb 2020 06:08:02 -0700
Subject: [PATCH 04/10] Add initial spec for trapz and trapz_weights
---
src/stdlib_experimental_quadrature.md | 79 +++++++++++++++++++++++++++
1 file changed, 79 insertions(+)
create mode 100644 src/stdlib_experimental_quadrature.md
diff --git a/src/stdlib_experimental_quadrature.md b/src/stdlib_experimental_quadrature.md
new file mode 100644
index 000000000..4627f2cb9
--- /dev/null
+++ b/src/stdlib_experimental_quadrature.md
@@ -0,0 +1,79 @@
+# 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` - compute 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 an 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 only 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
+
+```
From 0b385e745b192d87b08146e49a8f44ad7e8f0794 Mon Sep 17 00:00:00 2001
From: Nathaniel Shaffer <nrshaffer@protonmail.com>
Date: Wed, 5 Feb 2020 06:34:10 -0700
Subject: [PATCH 05/10] Add spec for simps
---
src/stdlib_experimental_quadrature.md | 36 +++++++++++++++++++++++++++
1 file changed, 36 insertions(+)
diff --git a/src/stdlib_experimental_quadrature.md b/src/stdlib_experimental_quadrature.md
index 4627f2cb9..fb4e3e6b8 100644
--- a/src/stdlib_experimental_quadrature.md
+++ b/src/stdlib_experimental_quadrature.md
@@ -77,3 +77,39 @@ program demo_trapz_weights
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 `u(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`.
+
+If the size of `y` is even, the result
+
+### Example
+
+TBD
From 770671d74413baae48bc830cc40a2872b7b7291c Mon Sep 17 00:00:00 2001
From: Nathaniel Shaffer <nrshaffer@protonmail.com>
Date: Thu, 6 Feb 2020 06:33:54 -0700
Subject: [PATCH 06/10] typos
---
src/stdlib_experimental_quadrature.md | 4 +---
1 file changed, 1 insertion(+), 3 deletions(-)
diff --git a/src/stdlib_experimental_quadrature.md b/src/stdlib_experimental_quadrature.md
index fb4e3e6b8..fc505d5f7 100644
--- a/src/stdlib_experimental_quadrature.md
+++ b/src/stdlib_experimental_quadrature.md
@@ -82,7 +82,7 @@ end program
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 `u(even:even+4)` and the ordinary Simpon's rule will be used elsewhere.
+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
@@ -108,8 +108,6 @@ 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 even, the result
-
### Example
TBD
From 8b6dcb7462f313f0f72844bc07e9d19440ddada2 Mon Sep 17 00:00:00 2001
From: Nathaniel Shaffer <nrshaffer@protonmail.com>
Date: Thu, 6 Feb 2020 06:42:23 -0700
Subject: [PATCH 07/10] minor wording changes
---
src/stdlib_experimental_quadrature.md | 6 +++---
1 file changed, 3 insertions(+), 3 deletions(-)
diff --git a/src/stdlib_experimental_quadrature.md b/src/stdlib_experimental_quadrature.md
index fc505d5f7..31f9b445f 100644
--- a/src/stdlib_experimental_quadrature.md
+++ b/src/stdlib_experimental_quadrature.md
@@ -44,7 +44,7 @@ program demo_trapz
end program
```
-## `trapz_weights` - compute trapezoidal rule weights for given abscissas
+## `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.
@@ -54,13 +54,13 @@ Given an array of abscissas `x`, computes the array of weights `w` such that if
### Arguments
-`x`: Shall be an array of type `real`.
+`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 only element of the result is zero.
+If the size of `x` is one, then the sole element of the result is zero.
### Example
From f7d732dc0315c17c86ba9494767bae32c1c521c4 Mon Sep 17 00:00:00 2001
From: Nathaniel Shaffer <nrshaffer@protonmail.com>
Date: Thu, 6 Feb 2020 06:42:31 -0700
Subject: [PATCH 08/10] add spec for simps_weights
---
src/stdlib_experimental_quadrature.md | 28 +++++++++++++++++++++++++++
1 file changed, 28 insertions(+)
diff --git a/src/stdlib_experimental_quadrature.md b/src/stdlib_experimental_quadrature.md
index 31f9b445f..c62638f64 100644
--- a/src/stdlib_experimental_quadrature.md
+++ b/src/stdlib_experimental_quadrature.md
@@ -111,3 +111,31 @@ If the size of `y` is two, the result is the same as if `trapz` had been called
### 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
From d1c91ecc56f244176d0fe386eb71dd6bdfdbae92 Mon Sep 17 00:00:00 2001
From: Nathaniel Shaffer <nrshaffer@protonmail.com>
Date: Fri, 7 Feb 2020 06:35:52 -0700
Subject: [PATCH 09/10] replace f.p. equality checks with tolerance checks
---
src/tests/quadrature/test_trapz.f90 | 84 ++++++++++++++---------------
1 file changed, 42 insertions(+), 42 deletions(-)
diff --git a/src/tests/quadrature/test_trapz.f90 b/src/tests/quadrature/test_trapz.f90
index b275643e8..7450878e0 100644
--- a/src/tests/quadrature/test_trapz.f90
+++ b/src/tests/quadrature/test_trapz.f90
@@ -12,7 +12,7 @@ program test_trapz
call test_trapz_weights_sp
call test_trapz_weights_dp
call test_trapz_weights_qp
-
+
call test_trapz_zero_sp
call test_trapz_zero_dp
call test_trapz_zero_qp
@@ -28,26 +28,26 @@ subroutine test_trapz_sp
integer :: i
print *, "test_trapz_sp"
-
+
y = [(real(i-1, sp), i = 1, n)]
val = trapz(y, 1.0_sp)
ans = 128.0_sp
- call assert(val == ans)
+ call assert(abs(val - ans) < epsilon(ans))
val = trapz(y, 0.5_sp)
ans = 64.0_sp
- call assert(val == ans)
+ call assert(abs(val - ans) < epsilon(ans))
x = [((i-1)*4.0_sp/real(n-1, sp), i = 1, n)]
val = trapz(y, x)
ans = 32.0_sp
- call assert(val == ans)
+ call assert(abs(val - ans) < epsilon(ans))
x = y**2
val = trapz(y, x)
ans = 2728.0_sp
- call assert(val == ans)
+ call assert(abs(val - ans) < epsilon(ans))
end subroutine test_trapz_sp
subroutine test_trapz_dp
@@ -59,26 +59,26 @@ subroutine test_trapz_dp
integer :: i
print *, "test_trapz_dp"
-
+
y = [(real(i-1, dp), i = 1, n)]
val = trapz(y, 1.0_dp)
ans = 128.0_dp
- call assert(val == ans)
+ call assert(abs(val - ans) < epsilon(ans))
val = trapz(y, 0.5_dp)
ans = 64.0_dp
- call assert(val == ans)
+ call assert(abs(val - ans) < epsilon(ans))
x = [((i-1)*4.0_dp/real(n-1, dp), i = 1, n)]
val = trapz(y, x)
ans = 32.0_dp
- call assert(val == ans)
+ call assert(abs(val - ans) < epsilon(ans))
x = y**2
val = trapz(y, x)
ans = 2728.0_sp
- call assert(val == ans)
+ call assert(abs(val - ans) < epsilon(ans))
end subroutine test_trapz_dp
@@ -91,26 +91,26 @@ subroutine test_trapz_qp
integer :: i
print *, "test_trapz_qp"
-
+
y = [(real(i-1, qp), i = 1, n)]
val = trapz(y, 1.0_qp)
ans = 128.0_qp
- call assert(val == ans)
+ call assert(abs(val - ans) < epsilon(ans))
val = trapz(y, 0.5_qp)
ans = 64.0_qp
- call assert(val == ans)
+ call assert(abs(val - ans) < epsilon(ans))
x = [((i-1)*4.0_qp/real(n-1, qp), i = 1, n)]
val = trapz(y, x)
ans = 32.0_qp
- call assert(val == ans)
+ call assert(abs(val - ans) < epsilon(ans))
x = y**2
val = trapz(y, x)
ans = 2728.0_qp
- call assert(val == ans)
+ call assert(abs(val - ans) < epsilon(ans))
end subroutine test_trapz_qp
@@ -124,20 +124,20 @@ subroutine test_trapz_weights_sp
real(sp) :: ans
print *, "test_trapz_weights_sp"
-
+
y = [(real(i-1, sp), i = 1, n)]
x = y
w = trapz_weights(x)
val = dot_product(w, y)
ans = trapz(y, x)
- call assert(val == ans)
+ call assert(abs(val - ans) < epsilon(ans))
x = y**2
w = trapz_weights(x)
val = dot_product(w, y)
ans = trapz(y, x)
- call assert(val == ans)
+ call assert(abs(val - ans) < epsilon(ans))
end subroutine test_trapz_weights_sp
@@ -152,24 +152,24 @@ subroutine test_trapz_weights_dp
real(dp) :: ans
print *, "test_trapz_weights_dp"
-
+
y = [(real(i-1, dp), i = 1, n)]
x = y
w = trapz_weights(x)
val = dot_product(w, y)
ans = trapz(y, x)
- call assert(val == ans)
+ call assert(abs(val - ans) < epsilon(ans))
x = y**2
w = trapz_weights(x)
val = dot_product(w, y)
ans = trapz(y, x)
- call assert(val == ans)
+ call assert(abs(val - ans) < epsilon(ans))
end subroutine test_trapz_weights_dp
-
+
subroutine test_trapz_weights_qp
integer, parameter :: n = 17
real(qp), dimension(n) :: y
@@ -180,20 +180,20 @@ subroutine test_trapz_weights_qp
real(qp) :: ans
print *, "test_trapz_weights_qp"
-
+
y = [(real(i-1, qp), i = 1, n)]
x = y
w = trapz_weights(x)
val = dot_product(w, y)
ans = trapz(y, x)
- call assert(val == ans)
+ call assert(abs(val - ans) < epsilon(ans))
x = y**2
w = trapz_weights(x)
val = dot_product(w, y)
ans = trapz(y, x)
- call assert(val == ans)
+ call assert(abs(val - ans) < epsilon(ans))
end subroutine test_trapz_weights_qp
@@ -202,11 +202,11 @@ subroutine test_trapz_zero_sp
real(sp), dimension(0) :: a
print *, "test_trapz_zero_sp"
-
- call assert(trapz(a, 1.0_sp) == 0.0_sp)
- call assert(trapz([1.0_sp], 1.0_sp) == 0.0_sp)
- call assert(trapz(a, a) == 0.0_sp)
- call assert(trapz([1.0_sp], [1.0_sp]) == 0.0_sp)
+
+ call assert(abs(trapz(a, 1.0_sp)) < epsilon(0.0_sp))
+ call assert(abs(trapz([1.0_sp], 1.0_sp)) < epsilon(0.0_sp))
+ call assert(abs(trapz(a, a)) < epsilon(0.0_sp))
+ call assert(abs(trapz([1.0_sp], [1.0_sp])) < epsilon(0.0_sp))
end subroutine test_trapz_zero_sp
@@ -214,11 +214,11 @@ subroutine test_trapz_zero_dp
real(dp), dimension(0) :: a
print *, "test_trapz_zero_dp"
-
- call assert(trapz(a, 1.0_dp) == 0.0_dp)
- call assert(trapz([1.0_dp], 1.0_dp) == 0.0_dp)
- call assert(trapz(a, a) == 0.0_dp)
- call assert(trapz([1.0_dp], [1.0_dp]) == 0.0_dp)
+
+ call assert(trapz(a, 1.0_dp)) < epsilon(0.0_dp))
+ call assert(abs(trapz([1.0_dp], 1.0_dp)) < epsilon(0.0_dp))
+ call assert(abs(trapz(a, a)) < epsilon(0.0_dp))
+ call assert(abs(trapz([1.0_dp], [1.0_dp])) < epsilon(0.0_dp))
end subroutine test_trapz_zero_dp
@@ -226,11 +226,11 @@ subroutine test_trapz_zero_qp
real(qp), dimension(0) :: a
print *, "test_trapz_zero_qp"
-
- call assert(trapz(a, 1.0_qp) == 0.0_qp)
- call assert(trapz([1.0_qp], 1.0_qp) == 0.0_qp)
- call assert(trapz(a, a) == 0.0_qp)
- call assert(trapz([1.0_qp], [1.0_qp]) == 0.0_qp)
+
+ call assert(trapz(a, 1.0_qp) < epsilon(0.0_qp))
+ call assert(trapz([1.0_qp], 1.0_qp) < epsilon(0.0_qp))
+ call assert(trapz(a, a) < epsilon(0.0_qp))
+ call assert(trapz([1.0_qp], [1.0_qp]) < epsilon(0.0_qp))
end subroutine test_trapz_zero_qp
-
+
end program test_trapz
From ccaf5485d1f07d5958a957083a6193a596aeb9bd Mon Sep 17 00:00:00 2001
From: Nathaniel Shaffer <nrshaffer@protonmail.com>
Date: Fri, 7 Feb 2020 06:41:31 -0700
Subject: [PATCH 10/10] fix compile error
---
src/tests/quadrature/test_trapz.f90 | 10 +++++-----
1 file changed, 5 insertions(+), 5 deletions(-)
diff --git a/src/tests/quadrature/test_trapz.f90 b/src/tests/quadrature/test_trapz.f90
index 7450878e0..1265e4b6b 100644
--- a/src/tests/quadrature/test_trapz.f90
+++ b/src/tests/quadrature/test_trapz.f90
@@ -215,7 +215,7 @@ subroutine test_trapz_zero_dp
print *, "test_trapz_zero_dp"
- call assert(trapz(a, 1.0_dp)) < epsilon(0.0_dp))
+ call assert(abs(trapz(a, 1.0_dp)) < epsilon(0.0_dp))
call assert(abs(trapz([1.0_dp], 1.0_dp)) < epsilon(0.0_dp))
call assert(abs(trapz(a, a)) < epsilon(0.0_dp))
call assert(abs(trapz([1.0_dp], [1.0_dp])) < epsilon(0.0_dp))
@@ -227,10 +227,10 @@ subroutine test_trapz_zero_qp
print *, "test_trapz_zero_qp"
- call assert(trapz(a, 1.0_qp) < epsilon(0.0_qp))
- call assert(trapz([1.0_qp], 1.0_qp) < epsilon(0.0_qp))
- call assert(trapz(a, a) < epsilon(0.0_qp))
- call assert(trapz([1.0_qp], [1.0_qp]) < epsilon(0.0_qp))
+ call assert(abs(trapz(a, 1.0_qp)) < epsilon(0.0_qp))
+ call assert(abs(trapz([1.0_qp], 1.0_qp)) < epsilon(0.0_qp))
+ call assert(abs(trapz(a, a)) < epsilon(0.0_qp))
+ call assert(abs(trapz([1.0_qp], [1.0_qp])) < epsilon(0.0_qp))
end subroutine test_trapz_zero_qp
end program test_trapz