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
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
|
subroutine dqk21(f,a,b,result,abserr,resabs,resasc)
c***begin prologue dqk21
c***date written 800101 (yymmdd)
c***revision date 830518 (yymmdd)
c***category no. h2a1a2
c***keywords 21-point gauss-kronrod rules
c***author piessens,robert,appl. math. & progr. div. - k.u.leuven
c de doncker,elise,appl. math. & progr. div. - k.u.leuven
c***purpose to compute i = integral of f over (a,b), with error
c estimate
c j = integral of abs(f) over (a,b)
c***description
c
c integration rules
c standard fortran subroutine
c double precision version
c
c parameters
c on entry
c f - double precision
c function subprogram defining the integrand
c function f(x). the actual name for f needs to be
c declared e x t e r n a l in the driver program.
c
c a - double precision
c lower limit of integration
c
c b - double precision
c upper limit of integration
c
c on return
c result - double precision
c approximation to the integral i
c result is computed by applying the 21-point
c kronrod rule (resk) obtained by optimal addition
c of abscissae to the 10-point gauss rule (resg).
c
c abserr - double precision
c estimate of the modulus of the absolute error,
c which should not exceed abs(i-result)
c
c resabs - double precision
c approximation to the integral j
c
c resasc - double precision
c approximation to the integral of abs(f-i/(b-a))
c over (a,b)
c
c***references (none)
c***routines called d1mach
c***end prologue dqk21
c
double precision a,absc,abserr,b,centr,dabs,dhlgth,dmax1,dmin1,
* d1mach,epmach,f,fc,fsum,fval1,fval2,fv1,fv2,hlgth,resabs,resasc,
* resg,resk,reskh,result,uflow,wg,wgk,xgk
integer j,jtw,jtwm1
external f
c
dimension fv1(10),fv2(10),wg(5),wgk(11),xgk(11)
c
c the abscissae and weights are given for the interval (-1,1).
c because of symmetry only the positive abscissae and their
c corresponding weights are given.
c
c xgk - abscissae of the 21-point kronrod rule
c xgk(2), xgk(4), ... abscissae of the 10-point
c gauss rule
c xgk(1), xgk(3), ... abscissae which are optimally
c added to the 10-point gauss rule
c
c wgk - weights of the 21-point kronrod rule
c
c wg - weights of the 10-point gauss rule
c
c
c gauss quadrature weights and kronron quadrature abscissae and weights
c as evaluated with 80 decimal digit arithmetic by l. w. fullerton,
c bell labs, nov. 1981.
c
data wg ( 1) / 0.0666713443 0868813759 3568809893 332 d0 /
data wg ( 2) / 0.1494513491 5058059314 5776339657 697 d0 /
data wg ( 3) / 0.2190863625 1598204399 5534934228 163 d0 /
data wg ( 4) / 0.2692667193 0999635509 1226921569 469 d0 /
data wg ( 5) / 0.2955242247 1475287017 3892994651 338 d0 /
c
data xgk ( 1) / 0.9956571630 2580808073 5527280689 003 d0 /
data xgk ( 2) / 0.9739065285 1717172007 7964012084 452 d0 /
data xgk ( 3) / 0.9301574913 5570822600 1207180059 508 d0 /
data xgk ( 4) / 0.8650633666 8898451073 2096688423 493 d0 /
data xgk ( 5) / 0.7808177265 8641689706 3717578345 042 d0 /
data xgk ( 6) / 0.6794095682 9902440623 4327365114 874 d0 /
data xgk ( 7) / 0.5627571346 6860468333 9000099272 694 d0 /
data xgk ( 8) / 0.4333953941 2924719079 9265943165 784 d0 /
data xgk ( 9) / 0.2943928627 0146019813 1126603103 866 d0 /
data xgk ( 10) / 0.1488743389 8163121088 4826001129 720 d0 /
data xgk ( 11) / 0.0000000000 0000000000 0000000000 000 d0 /
c
data wgk ( 1) / 0.0116946388 6737187427 8064396062 192 d0 /
data wgk ( 2) / 0.0325581623 0796472747 8818972459 390 d0 /
data wgk ( 3) / 0.0547558965 7435199603 1381300244 580 d0 /
data wgk ( 4) / 0.0750396748 1091995276 7043140916 190 d0 /
data wgk ( 5) / 0.0931254545 8369760553 5065465083 366 d0 /
data wgk ( 6) / 0.1093871588 0229764189 9210590325 805 d0 /
data wgk ( 7) / 0.1234919762 6206585107 7958109831 074 d0 /
data wgk ( 8) / 0.1347092173 1147332592 8054001771 707 d0 /
data wgk ( 9) / 0.1427759385 7706008079 7094273138 717 d0 /
data wgk ( 10) / 0.1477391049 0133849137 4841515972 068 d0 /
data wgk ( 11) / 0.1494455540 0291690566 4936468389 821 d0 /
c
c
c list of major variables
c -----------------------
c
c centr - mid point of the interval
c hlgth - half-length of the interval
c absc - abscissa
c fval* - function value
c resg - result of the 10-point gauss formula
c resk - result of the 21-point kronrod formula
c reskh - approximation to the mean value of f over (a,b),
c i.e. to i/(b-a)
c
c
c machine dependent constants
c ---------------------------
c
c epmach is the largest relative spacing.
c uflow is the smallest positive magnitude.
c
c***first executable statement dqk21
epmach = d1mach(4)
uflow = d1mach(1)
c
centr = 0.5d+00*(a+b)
hlgth = 0.5d+00*(b-a)
dhlgth = dabs(hlgth)
c
c compute the 21-point kronrod approximation to
c the integral, and estimate the absolute error.
c
resg = 0.0d+00
fc = f(centr)
resk = wgk(11)*fc
resabs = dabs(resk)
do 10 j=1,5
jtw = 2*j
absc = hlgth*xgk(jtw)
fval1 = f(centr-absc)
fval2 = f(centr+absc)
fv1(jtw) = fval1
fv2(jtw) = fval2
fsum = fval1+fval2
resg = resg+wg(j)*fsum
resk = resk+wgk(jtw)*fsum
resabs = resabs+wgk(jtw)*(dabs(fval1)+dabs(fval2))
10 continue
do 15 j = 1,5
jtwm1 = 2*j-1
absc = hlgth*xgk(jtwm1)
fval1 = f(centr-absc)
fval2 = f(centr+absc)
fv1(jtwm1) = fval1
fv2(jtwm1) = fval2
fsum = fval1+fval2
resk = resk+wgk(jtwm1)*fsum
resabs = resabs+wgk(jtwm1)*(dabs(fval1)+dabs(fval2))
15 continue
reskh = resk*0.5d+00
resasc = wgk(11)*dabs(fc-reskh)
do 20 j=1,10
resasc = resasc+wgk(j)*(dabs(fv1(j)-reskh)+dabs(fv2(j)-reskh))
20 continue
result = resk*hlgth
resabs = resabs*dhlgth
resasc = resasc*dhlgth
abserr = dabs((resk-resg)*hlgth)
if(resasc.ne.0.0d+00.and.abserr.ne.0.0d+00)
* abserr = resasc*dmin1(0.1d+01,(0.2d+03*abserr/resasc)**1.5d+00)
if(resabs.gt.uflow/(0.5d+02*epmach)) abserr = dmax1
* ((epmach*0.5d+02)*resabs,abserr)
return
end
|
