Skip to content

Commit 0cd354a

Browse files
jvdp1jvdp1
authored
Merge pull request #209 from nshaffer/dev-quadrature
Quadrature: spec tweaks & Simpson's rule
2 parents 4fbb208 + 9b6c3d0 commit 0cd354a

7 files changed

+879
-54
lines changed

‎doc/specs/stdlib_experimental_quadrature.md

+36-14
Original file line numberDiff line numberDiff line change
@@ -76,26 +76,26 @@ program demo_trapz_weights
7676
real :: x(5) = [0., 1., 2., 3., 4.]
7777
real :: y(5) = x**2
7878
real :: w(5)
79-
w = trapz_weight(x)
80-
print *, dot_product(w, y)
79+
w = trapz_weights(x)
80+
print *, sum(w*y)
8181
! 22.0
8282
end program demo_trapz_weights
8383
8484
```
8585

86-
## `simps` - integrate sampled values using Simpson's rule (to be implemented)
86+
## `simps` - integrate sampled values using Simpson's rule
8787

8888
### Description
8989

9090
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`.
9191

92-
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.
92+
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.
9393

9494
### Syntax
9595

96-
`result = simps(y, x [, even])`
96+
`result = [[stdlib_experimental_quadrature(module):simps(interface)]](y, x [, even])`
9797

98-
`result = simps(y, dx [, even])`
98+
`result = [[stdlib_experimental_quadrature(module):simps(interface)]](y, dx [, even])`
9999

100100
### Arguments
101101

@@ -105,37 +105,48 @@ Simpson's rule is defined for odd-length arrays only. For even-length arrays, an
105105

106106
`dx`: Shall be a scalar of type `real` having the same kind as `y`.
107107

108-
`even`: (Optional) Shall be a scalar integer of default kind. Its default value is `1`.
108+
`even`: (Optional) Shall be a default-kind `integer`.
109109

110110
### Return value
111111

112112
The result is a scalar of type `real` having the same kind as `y`.
113113

114114
If the size of `y` is zero or one, the result is zero.
115115

116-
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`.
116+
If the size of `y` is two, the result is the same as if `trapz` had been called instead.
117117

118118
### Example
119119

120-
TBD
120+
```fortran
121+
program demo_simps
122+
use stdlib_experimental_quadrature, only: simps
123+
implicit none
124+
real :: x(5) = [0., 1., 2., 3., 4.]
125+
real :: y(5) = 3.*x**2
126+
print *, simps(y, x)
127+
! 64.0
128+
print *, simps(y, 0.5)
129+
! 32.0
130+
end program demo_simps
131+
```
121132

122-
## `simps_weights` - Simpson's rule weights for given abscissas (to be implemented)
133+
## `simps_weights` - Simpson's rule weights for given abscissas
123134

124135
### Description
125136

126137
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.
127138

128-
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.
139+
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.
129140

130141
### Syntax
131142

132-
`result = simps_weights(x [, even])`
143+
`result = [[stdlib_experimental_quadrature(module):simps_weights(interface)]](x [, even])`
133144

134145
### Arguments
135146

136147
`x`: Shall be a rank-one array of type `real`.
137148

138-
`even`: (Optional) Shall be a scalar integer of default kind. Its default value is `1`.
149+
`even`: (Optional) Shall be a default-kind `integer`.
139150

140151
### Return value
141152

@@ -147,4 +158,15 @@ If the size of `x` is two, then the result is the same as if `trapz_weights` had
147158

148159
### Example
149160

150-
TBD
161+
```fortran
162+
program demo_simps_weights
163+
use stdlib_experimental_quadrature, only: simps_weights
164+
implicit none
165+
real :: x(5) = [0., 1., 2., 3., 4.]
166+
real :: y(5) = 3.*x**2
167+
real :: w(5)
168+
w = simps_weights(x)
169+
print *, sum(w*y)
170+
! 64.0
171+
end program demo_simps_weights
172+
```

‎src/CMakeLists.txt

+1
Original file line numberDiff line numberDiff line change
@@ -14,6 +14,7 @@ set(fppFiles
1414
stdlib_experimental_stats_var.fypp
1515
stdlib_experimental_quadrature.fypp
1616
stdlib_experimental_quadrature_trapz.fypp
17+
stdlib_experimental_quadrature_simps.fypp
1718
)
1819

1920

‎src/stdlib_experimental_quadrature.fypp

+53-35
Original file line numberDiff line numberDiff line change
@@ -1,5 +1,4 @@
11
#:include "common.fypp"
2-
32
module stdlib_experimental_quadrature
43
!! ([Specification](../page/specs/stdlib_experimental_quadrature.html#description))
54
use stdlib_experimental_kinds, only: sp, dp, qp
@@ -18,63 +17,82 @@ module stdlib_experimental_quadrature
1817
interface trapz
1918
!! Integrates sampled values using trapezoidal rule
2019
!! ([Specification](../page/specs/stdlib_experimental_quadrature.html#description))
21-
#:for KIND in REAL_KINDS
22-
pure module function trapz_dx_${KIND}$(y, dx) result(integral)
23-
real(${KIND}$), dimension(:), intent(in) :: y
24-
real(${KIND}$), intent(in) :: dx
25-
real(${KIND}$) :: integral
26-
end function trapz_dx_${KIND}$
20+
#:for k1, t1 in REAL_KINDS_TYPES
21+
pure module function trapz_dx_${k1}$(y, dx) result(integral)
22+
${t1}$, dimension(:), intent(in) :: y
23+
${t1}$, intent(in) :: dx
24+
${t1}$ :: integral
25+
end function trapz_dx_${k1}$
2726
#:endfor
28-
#:for KIND in REAL_KINDS
29-
pure module function trapz_x_${KIND}$(y, x) result(integral)
30-
real(${KIND}$), dimension(:), intent(in) :: y
31-
real(${KIND}$), dimension(size(y)), intent(in) :: x
32-
real(${KIND}$) :: integral
33-
end function trapz_x_${KIND}$
27+
#:for k1, t1 in REAL_KINDS_TYPES
28+
module function trapz_x_${k1}$(y, x) result(integral)
29+
${t1}$, dimension(:), intent(in) :: y
30+
${t1}$, dimension(:), intent(in) :: x
31+
${t1}$ :: integral
32+
end function trapz_x_${k1}$
3433
#:endfor
3534
end interface trapz
3635

3736

3837
interface trapz_weights
3938
!! Integrates sampled values using trapezoidal rule weights for given abscissas
4039
!! ([Specification](../page/specs/stdlib_experimental_quadrature.html#description_1))
41-
#:for KIND in REAL_KINDS
42-
pure module function trapz_weights_${KIND}$(x) result(w)
43-
real(${KIND}$), dimension(:), intent(in) :: x
44-
real(${KIND}$), dimension(size(x)) :: w
45-
end function trapz_weights_${KIND}$
40+
#:for k1, t1 in REAL_KINDS_TYPES
41+
pure module function trapz_weights_${k1}$(x) result(w)
42+
${t1}$, dimension(:), intent(in) :: x
43+
${t1}$, dimension(size(x)) :: w
44+
end function trapz_weights_${k1}$
4645
#:endfor
4746
end interface trapz_weights
4847

4948

5049
interface simps
51-
#:for KIND in REAL_KINDS
52-
pure module function simps_dx_${KIND}$(y, dx, even) result(integral)
53-
real(${KIND}$), dimension(:), intent(in) :: y
54-
real(${KIND}$), intent(in) :: dx
50+
!! Integrates sampled values using Simpson's rule
51+
!! ([Specification](../page/specs/stdlib_experimental_quadrature.html#description_3))
52+
! "recursive" is an implementation detail
53+
#:for k1, t1 in REAL_KINDS_TYPES
54+
pure recursive module function simps_dx_${k1}$(y, dx, even) result(integral)
55+
${t1}$, dimension(:), intent(in) :: y
56+
${t1}$, intent(in) :: dx
5557
integer, intent(in), optional :: even
56-
real(${KIND}$) :: integral
57-
end function simps_dx_${KIND}$
58+
${t1}$ :: integral
59+
end function simps_dx_${k1}$
5860
#:endfor
59-
#:for KIND in REAL_KINDS
60-
pure module function simps_x_${KIND}$(y, x, even) result(integral)
61-
real(${KIND}$), dimension(:), intent(in) :: y
62-
real(${KIND}$), dimension(size(y)), intent(in) :: x
61+
#:for k1, t1 in REAL_KINDS_TYPES
62+
recursive module function simps_x_${k1}$(y, x, even) result(integral)
63+
${t1}$, dimension(:), intent(in) :: y
64+
${t1}$, dimension(:), intent(in) :: x
6365
integer, intent(in), optional :: even
64-
real(${KIND}$) :: integral
65-
end function simps_x_${KIND}$
66+
${t1}$ :: integral
67+
end function simps_x_${k1}$
6668
#:endfor
6769
end interface simps
6870

6971

7072
interface simps_weights
71-
#:for KIND in REAL_KINDS
72-
pure module function simps_weights_${KIND}$(x, even) result(w)
73-
real(${KIND}$), dimension(:), intent(in) :: x
74-
real(${KIND}$), dimension(size(x)) :: w
73+
!! Integrates sampled values using trapezoidal rule weights for given abscissas
74+
!! ([Specification](../page/specs/stdlib_experimental_quadrature.html#description_3))
75+
#:for k1, t1 in REAL_KINDS_TYPES
76+
pure recursive module function simps_weights_${k1}$(x, even) result(w)
77+
${t1}$, dimension(:), intent(in) :: x
7578
integer, intent(in), optional :: even
76-
end function simps_weights_${KIND}$
79+
${t1}$, dimension(size(x)) :: w
80+
end function simps_weights_${k1}$
7781
#:endfor
7882
end interface simps_weights
7983

84+
85+
! Interface for a simple f(x)-style integrand function.
86+
! Could become fancier as we learn about the performance
87+
! ramifications of different ways to do callbacks.
88+
abstract interface
89+
#:for k1, t1 in REAL_KINDS_TYPES
90+
pure function integrand_${k1}$(x) result(f)
91+
import :: ${k1}$
92+
${t1}$, intent(in) :: x
93+
${t1}$ :: f
94+
end function integrand_${k1}$
95+
#:endfor
96+
end interface
97+
8098
end module stdlib_experimental_quadrature

0 commit comments

Comments
 (0)