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
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
|
subroutine dqk41(f,a,b,result,abserr,resabs,resasc)
c***begin prologue dqk41
c***date written 800101 (yymmdd)
c***revision date 830518 (yymmdd)
c***category no. h2a1a2
c***keywords 41-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 calling 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 41-point
c gauss-kronrod rule (resk) obtained by optimal
c addition of abscissae to the 20-point gauss
c 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 integal of abs(f-i/(b-a))
c over (a,b)
c
c***references (none)
c***routines called d1mach
c***end prologue dqk41
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(20),fv2(20),xgk(21),wgk(21),wg(10)
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 41-point gauss-kronrod rule
c xgk(2), xgk(4), ... abscissae of the 20-point
c gauss rule
c xgk(1), xgk(3), ... abscissae which are optimally
c added to the 20-point gauss rule
c
c wgk - weights of the 41-point gauss-kronrod rule
c
c wg - weights of the 20-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.0176140071 3915211831 1861962351 853 d0 /
data wg ( 2) / 0.0406014298 0038694133 1039952274 932 d0 /
data wg ( 3) / 0.0626720483 3410906356 9506535187 042 d0 /
data wg ( 4) / 0.0832767415 7670474872 4758143222 046 d0 /
data wg ( 5) / 0.1019301198 1724043503 6750135480 350 d0 /
data wg ( 6) / 0.1181945319 6151841731 2377377711 382 d0 /
data wg ( 7) / 0.1316886384 4917662689 8494499748 163 d0 /
data wg ( 8) / 0.1420961093 1838205132 9298325067 165 d0 /
data wg ( 9) / 0.1491729864 7260374678 7828737001 969 d0 /
data wg ( 10) / 0.1527533871 3072585069 8084331955 098 d0 /
c
data xgk ( 1) / 0.9988590315 8827766383 8315576545 863 d0 /
data xgk ( 2) / 0.9931285991 8509492478 6122388471 320 d0 /
data xgk ( 3) / 0.9815078774 5025025919 3342994720 217 d0 /
data xgk ( 4) / 0.9639719272 7791379126 7666131197 277 d0 /
data xgk ( 5) / 0.9408226338 3175475351 9982722212 443 d0 /
data xgk ( 6) / 0.9122344282 5132590586 7752441203 298 d0 /
data xgk ( 7) / 0.8782768112 5228197607 7442995113 078 d0 /
data xgk ( 8) / 0.8391169718 2221882339 4529061701 521 d0 /
data xgk ( 9) / 0.7950414288 3755119835 0638833272 788 d0 /
data xgk ( 10) / 0.7463319064 6015079261 4305070355 642 d0 /
data xgk ( 11) / 0.6932376563 3475138480 5490711845 932 d0 /
data xgk ( 12) / 0.6360536807 2651502545 2836696226 286 d0 /
data xgk ( 13) / 0.5751404468 1971031534 2946036586 425 d0 /
data xgk ( 14) / 0.5108670019 5082709800 4364050955 251 d0 /
data xgk ( 15) / 0.4435931752 3872510319 9992213492 640 d0 /
data xgk ( 16) / 0.3737060887 1541956067 2548177024 927 d0 /
data xgk ( 17) / 0.3016278681 1491300432 0555356858 592 d0 /
data xgk ( 18) / 0.2277858511 4164507808 0496195368 575 d0 /
data xgk ( 19) / 0.1526054652 4092267550 5220241022 678 d0 /
data xgk ( 20) / 0.0765265211 3349733375 4640409398 838 d0 /
data xgk ( 21) / 0.0000000000 0000000000 0000000000 000 d0 /
c
data wgk ( 1) / 0.0030735837 1852053150 1218293246 031 d0 /
data wgk ( 2) / 0.0086002698 5564294219 8661787950 102 d0 /
data wgk ( 3) / 0.0146261692 5697125298 3787960308 868 d0 /
data wgk ( 4) / 0.0203883734 6126652359 8010231432 755 d0 /
data wgk ( 5) / 0.0258821336 0495115883 4505067096 153 d0 /
data wgk ( 6) / 0.0312873067 7703279895 8543119323 801 d0 /
data wgk ( 7) / 0.0366001697 5820079803 0557240707 211 d0 /
data wgk ( 8) / 0.0416688733 2797368626 3788305936 895 d0 /
data wgk ( 9) / 0.0464348218 6749767472 0231880926 108 d0 /
data wgk ( 10) / 0.0509445739 2372869193 2707670050 345 d0 /
data wgk ( 11) / 0.0551951053 4828599474 4832372419 777 d0 /
data wgk ( 12) / 0.0591114008 8063957237 4967220648 594 d0 /
data wgk ( 13) / 0.0626532375 5478116802 5870122174 255 d0 /
data wgk ( 14) / 0.0658345971 3361842211 1563556969 398 d0 /
data wgk ( 15) / 0.0686486729 2852161934 5623411885 368 d0 /
data wgk ( 16) / 0.0710544235 5344406830 5790361723 210 d0 /
data wgk ( 17) / 0.0730306903 3278666749 5189417658 913 d0 /
data wgk ( 18) / 0.0745828754 0049918898 6581418362 488 d0 /
data wgk ( 19) / 0.0757044976 8455667465 9542775376 617 d0 /
data wgk ( 20) / 0.0763778676 7208073670 5502835038 061 d0 /
data wgk ( 21) / 0.0766007119 1799965644 5049901530 102 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 20-point gauss formula
c resk - result of the 41-point kronrod formula
c reskh - approximation to mean value of f over (a,b), i.e.
c to i/(b-a)
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 dqk41
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 41-point gauss-kronrod approximation to
c the integral, and estimate the absolute error.
c
resg = 0.0d+00
fc = f(centr)
resk = wgk(21)*fc
resabs = dabs(resk)
do 10 j=1,10
jtw = j*2
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,10
jtwm1 = j*2-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(21)*dabs(fc-reskh)
do 20 j=1,20
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.d+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
|
