1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
|
MODULE quadratures
IMPLICIT none
integer(kind=8) :: subintervals = 100
INTERFACE
FUNCTION integrate(ibeg, iend, myfun, p) result(val)
IMPLICIT none
! beginning of integration interval
real(kind=8), intent(in) :: ibeg
! ending of integration interval
real(kind=8), intent(in) :: iend
! function to be integrated
procedure(funint) :: myfun
! polynomial order
integer(kind=4), intent(in) :: p
! result of integration
real(kind=8) :: val
END FUNCTION integrate
END INTERFACE
INTERFACE
FUNCTION quadrature(qbeg, qend, fun) result(val)
IMPLICIT none
real(kind=8), intent(in) :: qbeg, qend
procedure(funint) :: fun
real(kind=8) :: val
END FUNCTION quadrature
END INTERFACE
INTERFACE
FUNCTION funint(x) result(y)
IMPLICIT none
real(kind=8), intent(in) :: x
real(kind=8) :: y
END FUNCTION funint
END INTERFACE
CONTAINS
FUNCTION newton_cotes(ibeg, iend, fun, p) result(val)
IMPLICIT none
real(kind=8), intent(in) :: ibeg
real(kind=8), intent(in) :: iend
procedure(funint) :: fun
integer(kind=4), intent(in) :: p
real(kind=8) :: val, subinterval_width, qbeg, qend
procedure(quadrature), pointer :: quad
integer(kind=8) :: i
subinterval_width = (iend - ibeg) / subintervals
SELECT CASE (p)
CASE (:-1)
STOP "negative interpolationg polynomial order passed"
CASE (0)
quad => rectangle
CASE (1)
quad => trapeze
CASE (2)
quad => simpson_1_third
CASE default
STOP "Newton-Cotes quadratures only implemented for order < 3"
END SELECT
! compute integral using quadrature pointed by quad
val = 0
DO i = 1, subintervals
qend = ibeg + i * subinterval_width
qbeg = ibeg + (i - 1) * subinterval_width
val = val + quad(qbeg, qend, fun)
END DO
END FUNCTION newton_cotes
FUNCTION rectangle(qbeg, qend, fun) result(val)
IMPLICIT none
real(kind=8), intent(in) :: qbeg, qend
procedure(funint) :: fun
real(kind=8) :: val
val = (qend - qbeg) * fun((qend + qbeg) * 0.5)
END FUNCTION rectangle
FUNCTION trapeze(qbeg, qend, fun) result(val)
IMPLICIT none
real(kind=8), intent(in) :: qbeg, qend
procedure(funint) :: fun
real(kind=8) :: val
val = (qend - qbeg) * 0.5 * (fun(qbeg) + fun(qend))
END FUNCTION trapeze
FUNCTION simpson_1_third(qbeg, qend, fun) result(val)
IMPLICIT none
real(kind=8), intent(in) :: qbeg, qend
procedure(funint) :: fun
real(kind=8) :: val
val = (qend - qbeg) * (1/6.0) * &
(fun(qbeg) + 4 * fun ((qbeg + qend) * 0.5) + fun(qend))
END FUNCTION simpson_1_third
END MODULE quadratures
|
