Skip to content

Commit cd4ef20

Browse files
nshaffernshaffer
committed
implement Simpson's rule
1 parent d1db8a7 commit cd4ef20

File tree

4 files changed

+784
-1
lines changed

4 files changed

+784
-1
lines changed

‎src/CMakeLists.txt

Lines changed: 1 addition & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -7,6 +7,7 @@ set(fppFiles
77
stdlib_experimental_stats_mean.fypp
88
stdlib_experimental_quadrature.fypp
99
stdlib_experimental_quadrature_trapz.fypp
10+
stdlib_experimental_quadrature_simps.fypp
1011
)
1112

1213

‎src/stdlib_experimental_quadrature_simps.fypp

Lines changed: 277 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,277 @@
1+
#:include "common.fypp"
2+
3+
submodule (stdlib_experimental_quadrature) stdlib_experimental_quadrature_simps
4+
5+
implicit none
6+
7+
! internal use only
8+
interface simps38
9+
#:for k1, t1 in REAL_KINDS_TYPES
10+
module procedure simps38_dx_${k1}$
11+
module procedure simps38_x_${k1}$
12+
#:endfor
13+
end interface simps38
14+
15+
! internal use only
16+
interface simps38_weights
17+
#:for k1, t1 in REAL_KINDS_TYPES
18+
module procedure simps38_weights_${k1}$
19+
#:endfor
20+
end interface simps38_weights
21+
22+
contains
23+
24+
#:for k1, t1 in REAL_KINDS_TYPES
25+
26+
pure recursive module function simps_dx_${k1}$(y, dx, even) result(integral)
27+
${t1}$, dimension(:), intent(in) :: y
28+
${t1}$, intent(in) :: dx
29+
integer, intent(in), optional :: even
30+
${t1}$ :: integral
31+
32+
integer :: n
33+
34+
n = size(y)
35+
36+
select case (n)
37+
case (0:1)
38+
integral = 0.0_${k1}$
39+
case (2)
40+
integral = 0.5_${k1}$*dx*(y(1) + y(2))
41+
case (3)
42+
integral = dx/3.0_${k1}$*(y(1) + 4*y(2) + y(3))
43+
case (4)
44+
integral = simps38(y, dx)
45+
! case (5) not needed; handled by default
46+
case (6) ! needs special handling because of averaged 3/8's rule case
47+
if (present(even)) then
48+
if (even < 0) then
49+
! 3/8 rule on left
50+
integral = simps38(y(1:4), dx) + simps(y(4:6), dx)
51+
return
52+
else if (even > 0) then
53+
! 3/8 rule on right
54+
integral = simps(y(1:3), dx) + simps38(y(3:6), dx)
55+
return
56+
else
57+
! fall through
58+
end if
59+
end if
60+
! either `even` not present or is zero
61+
! equivalent to averaging left and right
62+
integral = dx/48.0_${k1}$ * (17*(y(1) + y(6)) + 59*(y(2) + y(5)) + 44*(y(3) + y(4)))
63+
case default
64+
if (mod(n, 2) == 1) then
65+
integral = dx/3.0_${k1}$*(y(1) + 4*sum(y(2:n-1:2)) + 2*sum(y(3:n-2:2)) + y(n))
66+
else
67+
if (present(even)) then
68+
if (even < 0) then
69+
! 3/8th rule on left
70+
integral = simps38(y(1:4), dx) + simps(y(4:n), dx)
71+
return
72+
else if (even > 0) then
73+
! 3/8 rule on right
74+
integral = simps(y(1:n-3), dx) + simps38(y(n-3:n), dx)
75+
return
76+
else
77+
! fall through
78+
end if
79+
end if
80+
! either `even` not present or is zero
81+
! equivalent to averaging left and right
82+
integral = dx/48.0_${k1}$ * (17*(y(1) + y(n)) + 59*(y(2) + y(n-1)) &
83+
+ 43*(y(3) + y(n-2)) + 49*(y(4) + y(n-3)) + 48*sum(y(5:n-4)))
84+
end if
85+
end select
86+
end function simps_dx_${k1}$
87+
88+
#:endfor
89+
#:for k1, t1 in REAL_KINDS_TYPES
90+
91+
pure recursive module function simps_x_${k1}$(y, x, even) result(integral)
92+
${t1}$, dimension(:), intent(in) :: y
93+
${t1}$, dimension(size(y)), intent(in) :: x
94+
integer, intent(in), optional :: even
95+
${t1}$ :: integral
96+
97+
integer :: i
98+
integer :: n
99+
100+
${t1}$ :: h1, h2, h3
101+
${t1}$ :: a, b, c, d
102+
103+
n = size(y)
104+
105+
select case (n)
106+
case (0:1)
107+
integral = 0.0_${k1}$
108+
case (2)
109+
integral = 0.5_${k1}$*(x(2) - x(1))*(y(1) + y(2))
110+
case (3)
111+
h1 = x(2) - x(1)
112+
h2 = x(3) - x(2)
113+
a = (2*h1**2 + h1*h2 - h2**2)/(6*h1)
114+
b = (h1+h2)**3/(6*h1*h2)
115+
c = (2*h2**2 + h1*h2 - h1**2)/(6*h2)
116+
integral = a*y(1) + b*y(2) + c*y(3)
117+
case (4)
118+
integral = simps38(y, x)
119+
! case (6) unneeded?
120+
case default
121+
if (mod(n, 2) == 1) then
122+
integral = 0.0_${k1}$
123+
do i = 1, n-2, 2
124+
h1 = x(i+1) - x(i)
125+
h2 = x(i+2) - x(i+1)
126+
a = (2*h1**2 + h1*h2 - h2**2)/(6*h1)
127+
b = (h1+h2)**3/(6*h1*h2)
128+
c = (2*h2**2 + h1*h2 - h1**2)/(6*h2)
129+
integral = integral + a*y(i) + b*y(i+1) + c*y(i+2)
130+
end do
131+
else
132+
if (present(even)) then
133+
if (even < 0) then
134+
! 3/8 rule on left
135+
integral = simps38(y(1:4), x(1:4)) + simps(y(4:n), x(4:n))
136+
return
137+
else if (even > 0) then
138+
! 3/8 rule on right
139+
integral = simps(y(1:n-3), x(1:n-3)) + simps38(y(n-3:n), x(n-3:n))
140+
return
141+
else
142+
! fall through
143+
end if
144+
end if
145+
! either `even` not present or is zero
146+
integral = 0.5_${k1}$ * ( simps38(y(1:4), x(1:4)) + simps(y(4:n), x(4:n)) &
147+
+ simps(y(1:n-3), x(1:n-3)) + simps38(y(n-3:n), x(n-3:n)) )
148+
end if
149+
end select
150+
end function simps_x_${k1}$
151+
152+
#:endfor
153+
#:for k1, t1 in REAL_KINDS_TYPES
154+
155+
pure recursive module function simps_weights_${k1}$(x, even) result(w)
156+
${t1}$, dimension(:), intent(in) :: x
157+
integer, intent(in), optional :: even
158+
${t1}$, dimension(size(x)) :: w
159+
160+
integer :: i, n
161+
${t1}$ :: h1, h2, h3
162+
163+
n = size(x)
164+
165+
select case (n)
166+
case (0)
167+
! no action needed
168+
case (1)
169+
w(1) = 0.0_${k1}$
170+
case (2)
171+
w = 0.5_${k1}$*(x(2) - x(1))
172+
case (3)
173+
h1 = x(2) - x(1)
174+
h2 = x(3) - x(2)
175+
w(1) = (2*h1**2 + h1*h2 - h2**2)/(6*h1)
176+
w(2) = (h1+h2)**3/(6*h1*h2)
177+
w(3) = (2*h2**2 + h1*h2 - h1**2)/(6*h2)
178+
case (4)
179+
w = simps38_weights(x)
180+
case default
181+
if (mod(n, 2) == 1) then
182+
w = 0.0_${k1}$
183+
do i = 1, n-2, 2
184+
h1 = x(i+1) - x(i)
185+
h2 = x(i+2) - x(i+1)
186+
w(i) = w(i) + (2*h1**2 + h1*h2 - h2**2)/(6*h1)
187+
w(i+1) = w(i+1) + (h1+h2)**3/(6*h1*h2)
188+
w(i+2) = w(i+2) + (2*h2**2 + h1*h2 - h1**2)/(6*h2)
189+
end do
190+
else
191+
if (present(even)) then
192+
if (even < 0) then
193+
! 3/8 rule on left
194+
w = 0.0_${k1}$
195+
w(1:4) = simps38_weights(x(1:4))
196+
w(4:n) = w(4:n) + simps_weights(x(4:n)) ! position 4 needs both rules
197+
return
198+
else if (even > 0) then
199+
! 3/8 rule on right
200+
w = 0.0_${k1}$
201+
w(1:n-3) = simps_weights(x(1:n-3))
202+
w(n-3:n) = w(n-3:n) + simps38_weights(x(n-3:n)) ! position n-3 needs both rules
203+
return
204+
else
205+
! fall through
206+
end if
207+
end if
208+
! either `even` not present or is zero
209+
w = 0.0_${k1}$
210+
! 3/8 rule on left
211+
w(1:4) = simps38_weights(x(1:4))
212+
w(4:n) = w(4:n) + simps_weights(x(4:n))
213+
! 3/8 rule on right
214+
w(1:n-3) = w(1:n-3) + simps_weights(x(1:n-3))
215+
w(n-3:n) = w(n-3:n) + simps38_weights(x(n-3:n))
216+
! average
217+
w = 0.5_${k1}$ * w
218+
end if
219+
end select
220+
end function simps_weights_${k1}$
221+
222+
#:endfor
223+
#:for k1, t1 in REAL_KINDS_TYPES
224+
225+
pure function simps38_dx_${k1}$(y, dx) result (integral)
226+
${t1}$, dimension(4), intent(in) :: y
227+
${t1}$, intent(in) :: dx
228+
${t1}$ :: integral
229+
230+
integral = 3.0_${k1}$*dx/8.0_${k1}$ * (y(1) + y(4) + 3*(y(2) + y(3)))
231+
end function simps38_dx_${k1}$
232+
233+
#:endfor
234+
#: for k1, t1 in REAL_KINDS_TYPES
235+
236+
pure function simps38_x_${k1}$(y, x) result(integral)
237+
${t1}$, dimension(4), intent(in) :: y
238+
${t1}$, dimension(4), intent(in) :: x
239+
${t1}$ :: integral
240+
241+
${t1}$ :: h1, h2, h3
242+
${t1}$ :: a, b, c, d
243+
244+
h1 = x(2) - x(1)
245+
h2 = x(3) - x(2)
246+
h3 = x(4) - x(3)
247+
248+
a = (h1+h2+h3)*(3*h1**2 + 2*h1*h2 - 2*h1*h3 - h2**2 + h3**2)/(12*h1*(h1+h2))
249+
b = (h1+h2-h3)*(h1+h2+h3)**3/(12*h1*h2*(h2+h3))
250+
c = (h2+h3-h1)*(h1+h2+h3)**3/(12*h2*h3*(h1+h2))
251+
d = (h1+h2+h3)*(3*h3**2 + 2*h2*h3 - 2*h1*h3 - h2**2 + h1**2)/(12*h3*(h2+h3))
252+
253+
integral = a*y(1) + b*y(2) + c*y(3) + d*y(4)
254+
end function simps38_x_${k1}$
255+
256+
#:endfor
257+
#:for k1, t1 in REAL_KINDS_TYPES
258+
259+
pure function simps38_weights_${k1}$(x) result(w)
260+
${t1}$, intent(in) :: x(4)
261+
${t1}$ :: w(size(x))
262+
263+
${t1}$ :: h1, h2, h3
264+
265+
h1 = x(2) - x(1)
266+
h2 = x(3) - x(2)
267+
h3 = x(4) - x(3)
268+
269+
w(1) = (h1+h2+h3)*(3*h1**2 + 2*h1*h2 - 2*h1*h3 - h2**2 + h3**2)/(12*h1*(h1+h2))
270+
w(2) = (h1+h2-h3)*(h1+h2+h3)**3/(12*h1*h2*(h2+h3))
271+
w(3) = (h2+h3-h1)*(h1+h2+h3)**3/(12*h2*h3*(h1+h2))
272+
w(4) = (h1+h2+h3)*(3*h3**2 + 2*h2*h3 - 2*h1*h3 - h2**2 + h1**2)/(12*h3*(h2+h3))
273+
end function simps38_weights_${k1}$
274+
275+
#:endfor
276+
277+
end submodule stdlib_experimental_quadrature_simps

‎src/tests/quadrature/CMakeLists.txt

Lines changed: 1 addition & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -1,2 +1,2 @@
11
ADDTEST(trapz)
2-
2+
ADDTEST(simps)

0 commit comments

Comments
 (0)