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
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
|
subroutine dqng(f,a,b,epsabs,epsrel,result,abserr,neval,ier)
c***begin prologue dqng
c***date written 800101 (yymmdd)
c***revision date 810101 (yymmdd)
c***category no. h2a1a1
c***keywords automatic integrator, smooth integrand,
c non-adaptive, gauss-kronrod(patterson)
c***author piessens,robert,appl. math. & progr. div. - k.u.leuven
c de doncker,elise,appl math & progr. div. - k.u.leuven
c kahaner,david,nbs - modified (2/82)
c***purpose the routine calculates an approximation result to a
c given definite integral i = integral of f over (a,b),
c hopefully satisfying following claim for accuracy
c abs(i-result).le.max(epsabs,epsrel*abs(i)).
c***description
c
c non-adaptive integration
c standard fortran subroutine
c double precision version
c
c f - double precision
c function subprogram defining the integrand function
c f(x). the actual name for f needs to be declared
c 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 epsabs - double precision
c absolute accuracy requested
c epsrel - double precision
c relative accuracy requested
c if epsabs.le.0
c and epsrel.lt.max(50*rel.mach.acc.,0.5d-28),
c the routine will end with ier = 6.
c
c on return
c result - double precision
c approximation to the integral i
c result is obtained by applying the 21-point
c gauss-kronrod rule (res21) obtained by optimal
c addition of abscissae to the 10-point gauss rule
c (res10), or by applying the 43-point rule (res43)
c obtained by optimal addition of abscissae to the
c 21-point gauss-kronrod rule, or by applying the
c 87-point rule (res87) obtained by optimal addition
c of abscissae to the 43-point rule.
c
c abserr - double precision
c estimate of the modulus of the absolute error,
c which should equal or exceed abs(i-result)
c
c neval - integer
c number of integrand evaluations
c
c ier - ier = 0 normal and reliable termination of the
c routine. it is assumed that the requested
c accuracy has been achieved.
c ier.gt.0 abnormal termination of the routine. it is
c assumed that the requested accuracy has
c not been achieved.
c error messages
c ier = 1 the maximum number of steps has been
c executed. the integral is probably too
c difficult to be calculated by dqng.
c = 6 the input is invalid, because
c epsabs.le.0 and
c epsrel.lt.max(50*rel.mach.acc.,0.5d-28).
c result, abserr and neval are set to zero.
c
c***references (none)
c***routines called d1mach,xerror
c***end prologue dqng
c
double precision a,absc,abserr,b,centr,dabs,dhlgth,dmax1,dmin1,
* d1mach,epmach,epsabs,epsrel,f,fcentr,fval,fval1,fval2,fv1,fv2,
* fv3,fv4,hlgth,result,res10,res21,res43,res87,resabs,resasc,
* reskh,savfun,uflow,w10,w21a,w21b,w43a,w43b,w87a,w87b,x1,x2,x3,x4
integer ier,ipx,k,l,neval
external f
c
dimension fv1(5),fv2(5),fv3(5),fv4(5),x1(5),x2(5),x3(11),x4(22),
* w10(5),w21a(5),w21b(6),w43a(10),w43b(12),w87a(21),w87b(23),
* savfun(21)
c
c the following data statements contain the
c abscissae and weights of the integration rules used.
c
c x1 abscissae common to the 10-, 21-, 43- and 87-
c point rule
c x2 abscissae common to the 21-, 43- and 87-point rule
c x3 abscissae common to the 43- and 87-point rule
c x4 abscissae of the 87-point rule
c w10 weights of the 10-point formula
c w21a weights of the 21-point formula for abscissae x1
c w21b weights of the 21-point formula for abscissae x2
c w43a weights of the 43-point formula for abscissae x1, x3
c w43b weights of the 43-point formula for abscissae x3
c w87a weights of the 87-point formula for abscissae x1,
c x2, x3
c w87b weights of the 87-point formula for abscissae x4
c
c
c gauss-kronrod-patterson quadrature coefficients for use in
c quadpack routine qng. these coefficients were calculated with
c 101 decimal digit arithmetic by l. w. fullerton, bell labs, nov 1981.
c
data x1 ( 1) / 0.9739065285 1717172007 7964012084 452 d0 /
data x1 ( 2) / 0.8650633666 8898451073 2096688423 493 d0 /
data x1 ( 3) / 0.6794095682 9902440623 4327365114 874 d0 /
data x1 ( 4) / 0.4333953941 2924719079 9265943165 784 d0 /
data x1 ( 5) / 0.1488743389 8163121088 4826001129 720 d0 /
data w10 ( 1) / 0.0666713443 0868813759 3568809893 332 d0 /
data w10 ( 2) / 0.1494513491 5058059314 5776339657 697 d0 /
data w10 ( 3) / 0.2190863625 1598204399 5534934228 163 d0 /
data w10 ( 4) / 0.2692667193 0999635509 1226921569 469 d0 /
data w10 ( 5) / 0.2955242247 1475287017 3892994651 338 d0 /
c
data x2 ( 1) / 0.9956571630 2580808073 5527280689 003 d0 /
data x2 ( 2) / 0.9301574913 5570822600 1207180059 508 d0 /
data x2 ( 3) / 0.7808177265 8641689706 3717578345 042 d0 /
data x2 ( 4) / 0.5627571346 6860468333 9000099272 694 d0 /
data x2 ( 5) / 0.2943928627 0146019813 1126603103 866 d0 /
data w21a ( 1) / 0.0325581623 0796472747 8818972459 390 d0 /
data w21a ( 2) / 0.0750396748 1091995276 7043140916 190 d0 /
data w21a ( 3) / 0.1093871588 0229764189 9210590325 805 d0 /
data w21a ( 4) / 0.1347092173 1147332592 8054001771 707 d0 /
data w21a ( 5) / 0.1477391049 0133849137 4841515972 068 d0 /
data w21b ( 1) / 0.0116946388 6737187427 8064396062 192 d0 /
data w21b ( 2) / 0.0547558965 7435199603 1381300244 580 d0 /
data w21b ( 3) / 0.0931254545 8369760553 5065465083 366 d0 /
data w21b ( 4) / 0.1234919762 6206585107 7958109831 074 d0 /
data w21b ( 5) / 0.1427759385 7706008079 7094273138 717 d0 /
data w21b ( 6) / 0.1494455540 0291690566 4936468389 821 d0 /
c
data x3 ( 1) / 0.9993333609 0193208139 4099323919 911 d0 /
data x3 ( 2) / 0.9874334029 0808886979 5961478381 209 d0 /
data x3 ( 3) / 0.9548079348 1426629925 7919200290 473 d0 /
data x3 ( 4) / 0.9001486957 4832829362 5099494069 092 d0 /
data x3 ( 5) / 0.8251983149 8311415084 7066732588 520 d0 /
data x3 ( 6) / 0.7321483889 8930498261 2354848755 461 d0 /
data x3 ( 7) / 0.6228479705 3772523864 1159120344 323 d0 /
data x3 ( 8) / 0.4994795740 7105649995 2214885499 755 d0 /
data x3 ( 9) / 0.3649016613 4658076804 3989548502 644 d0 /
data x3 ( 10) / 0.2222549197 7660129649 8260928066 212 d0 /
data x3 ( 11) / 0.0746506174 6138332204 3914435796 506 d0 /
data w43a ( 1) / 0.0162967342 8966656492 4281974617 663 d0 /
data w43a ( 2) / 0.0375228761 2086950146 1613795898 115 d0 /
data w43a ( 3) / 0.0546949020 5825544214 7212685465 005 d0 /
data w43a ( 4) / 0.0673554146 0947808607 5553166302 174 d0 /
data w43a ( 5) / 0.0738701996 3239395343 2140695251 367 d0 /
data w43a ( 6) / 0.0057685560 5976979618 4184327908 655 d0 /
data w43a ( 7) / 0.0273718905 9324884208 1276069289 151 d0 /
data w43a ( 8) / 0.0465608269 1042883074 3339154433 824 d0 /
data w43a ( 9) / 0.0617449952 0144256449 6240336030 883 d0 /
data w43a ( 10) / 0.0713872672 6869339776 8559114425 516 d0 /
data w43b ( 1) / 0.0018444776 4021241410 0389106552 965 d0 /
data w43b ( 2) / 0.0107986895 8589165174 0465406741 293 d0 /
data w43b ( 3) / 0.0218953638 6779542810 2523123075 149 d0 /
data w43b ( 4) / 0.0325974639 7534568944 3882222526 137 d0 /
data w43b ( 5) / 0.0421631379 3519181184 7627924327 955 d0 /
data w43b ( 6) / 0.0507419396 0018457778 0189020092 084 d0 /
data w43b ( 7) / 0.0583793955 4261924837 5475369330 206 d0 /
data w43b ( 8) / 0.0647464049 5144588554 4689259517 511 d0 /
data w43b ( 9) / 0.0695661979 1235648452 8633315038 405 d0 /
data w43b ( 10) / 0.0728244414 7183320815 0939535192 842 d0 /
data w43b ( 11) / 0.0745077510 1417511827 3571813842 889 d0 /
data w43b ( 12) / 0.0747221475 1740300559 4425168280 423 d0 /
c
data x4 ( 1) / 0.9999029772 6272923449 0529830591 582 d0 /
data x4 ( 2) / 0.9979898959 8667874542 7496322365 960 d0 /
data x4 ( 3) / 0.9921754978 6068722280 8523352251 425 d0 /
data x4 ( 4) / 0.9813581635 7271277357 1916941623 894 d0 /
data x4 ( 5) / 0.9650576238 5838461912 8284110607 926 d0 /
data x4 ( 6) / 0.9431676131 3367059681 6416634507 426 d0 /
data x4 ( 7) / 0.9158064146 8550720959 1826430720 050 d0 /
data x4 ( 8) / 0.8832216577 7131650137 2117548744 163 d0 /
data x4 ( 9) / 0.8457107484 6241566660 5902011504 855 d0 /
data x4 ( 10) / 0.8035576580 3523098278 8739474980 964 d0 /
data x4 ( 11) / 0.7570057306 8549555832 8942793432 020 d0 /
data x4 ( 12) / 0.7062732097 8732181982 4094274740 840 d0 /
data x4 ( 13) / 0.6515894665 0117792253 4422205016 736 d0 /
data x4 ( 14) / 0.5932233740 5796108887 5273770349 144 d0 /
data x4 ( 15) / 0.5314936059 7083193228 5268948562 671 d0 /
data x4 ( 16) / 0.4667636230 4202284487 1966781659 270 d0 /
data x4 ( 17) / 0.3994248478 5921880473 2101665817 923 d0 /
data x4 ( 18) / 0.3298748771 0618828826 5053371824 597 d0 /
data x4 ( 19) / 0.2585035592 0216155180 2280975429 025 d0 /
data x4 ( 20) / 0.1856953965 6834665201 5917141167 606 d0 /
data x4 ( 21) / 0.1118422131 7990746817 2398359241 362 d0 /
data x4 ( 22) / 0.0373521233 9461987081 4998165437 704 d0 /
data w87a ( 1) / 0.0081483773 8414917290 0002878448 190 d0 /
data w87a ( 2) / 0.0187614382 0156282224 3935059003 794 d0 /
data w87a ( 3) / 0.0273474510 5005228616 1582829741 283 d0 /
data w87a ( 4) / 0.0336777073 1163793004 6581056957 588 d0 /
data w87a ( 5) / 0.0369350998 2042790761 4589586742 499 d0 /
data w87a ( 6) / 0.0028848724 3021153050 1334156248 695 d0 /
data w87a ( 7) / 0.0136859460 2271270188 8950035273 128 d0 /
data w87a ( 8) / 0.0232804135 0288831112 3409291030 404 d0 /
data w87a ( 9) / 0.0308724976 1171335867 5466394126 442 d0 /
data w87a ( 10) / 0.0356936336 3941877071 9351355457 044 d0 /
data w87a ( 11) / 0.0009152833 4520224136 0843392549 948 d0 /
data w87a ( 12) / 0.0053992802 1930047136 7738743391 053 d0 /
data w87a ( 13) / 0.0109476796 0111893113 4327826856 808 d0 /
data w87a ( 14) / 0.0162987316 9678733526 2665703223 280 d0 /
data w87a ( 15) / 0.0210815688 8920383511 2433060188 190 d0 /
data w87a ( 16) / 0.0253709697 6925382724 3467999831 710 d0 /
data w87a ( 17) / 0.0291896977 5647575250 1446154084 920 d0 /
data w87a ( 18) / 0.0323732024 6720278968 5788194889 595 d0 /
data w87a ( 19) / 0.0347830989 5036514275 0781997949 596 d0 /
data w87a ( 20) / 0.0364122207 3135178756 2801163687 577 d0 /
data w87a ( 21) / 0.0372538755 0304770853 9592001191 226 d0 /
data w87b ( 1) / 0.0002741455 6376207235 0016527092 881 d0 /
data w87b ( 2) / 0.0018071241 5505794294 8341311753 254 d0 /
data w87b ( 3) / 0.0040968692 8275916486 4458070683 480 d0 /
data w87b ( 4) / 0.0067582900 5184737869 9816577897 424 d0 /
data w87b ( 5) / 0.0095499576 7220164653 6053581325 377 d0 /
data w87b ( 6) / 0.0123294476 5224485369 4626639963 780 d0 /
data w87b ( 7) / 0.0150104473 4638895237 6697286041 943 d0 /
data w87b ( 8) / 0.0175489679 8624319109 9665352925 900 d0 /
data w87b ( 9) / 0.0199380377 8644088820 2278192730 714 d0 /
data w87b ( 10) / 0.0221949359 6101228679 6332102959 499 d0 /
data w87b ( 11) / 0.0243391471 2600080547 0360647041 454 d0 /
data w87b ( 12) / 0.0263745054 1483920724 1503786552 615 d0 /
data w87b ( 13) / 0.0282869107 8877120065 9968002987 960 d0 /
data w87b ( 14) / 0.0300525811 2809269532 2521110347 341 d0 /
data w87b ( 15) / 0.0316467513 7143992940 4586051078 883 d0 /
data w87b ( 16) / 0.0330504134 1997850329 0785944862 689 d0 /
data w87b ( 17) / 0.0342550997 0422606178 7082821046 821 d0 /
data w87b ( 18) / 0.0352624126 6015668103 3782717998 428 d0 /
data w87b ( 19) / 0.0360769896 2288870118 5500318003 895 d0 /
data w87b ( 20) / 0.0366986044 9845609449 8018047441 094 d0 /
data w87b ( 21) / 0.0371205492 6983257611 4119958413 599 d0 /
data w87b ( 22) / 0.0373342287 5193504032 1235449094 698 d0 /
data w87b ( 23) / 0.0373610737 6267902341 0321241766 599 d0 /
c
c list of major variables
c -----------------------
c
c centr - mid point of the integration interval
c hlgth - half-length of the integration interval
c fcentr - function value at mid point
c absc - abscissa
c fval - function value
c savfun - array of function values which have already been
c computed
c res10 - 10-point gauss result
c res21 - 21-point kronrod result
c res43 - 43-point result
c res87 - 87-point result
c resabs - approximation to the integral of abs(f)
c resasc - approximation to the integral of abs(f-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 dqng
epmach = d1mach(4)
uflow = d1mach(1)
c
c test on validity of parameters
c ------------------------------
c
result = 0.0d+00
abserr = 0.0d+00
neval = 0
ier = 6
if(epsabs.le.0.0d+00.and.epsrel.lt.dmax1(0.5d+02*epmach,0.5d-28))
* go to 80
hlgth = 0.5d+00*(b-a)
dhlgth = dabs(hlgth)
centr = 0.5d+00*(b+a)
fcentr = f(centr)
neval = 21
ier = 1
c
c compute the integral using the 10- and 21-point formula.
c
do 70 l = 1,3
go to (5,25,45),l
5 res10 = 0.0d+00
res21 = w21b(6)*fcentr
resabs = w21b(6)*dabs(fcentr)
do 10 k=1,5
absc = hlgth*x1(k)
fval1 = f(centr+absc)
fval2 = f(centr-absc)
fval = fval1+fval2
res10 = res10+w10(k)*fval
res21 = res21+w21a(k)*fval
resabs = resabs+w21a(k)*(dabs(fval1)+dabs(fval2))
savfun(k) = fval
fv1(k) = fval1
fv2(k) = fval2
10 continue
ipx = 5
do 15 k=1,5
ipx = ipx+1
absc = hlgth*x2(k)
fval1 = f(centr+absc)
fval2 = f(centr-absc)
fval = fval1+fval2
res21 = res21+w21b(k)*fval
resabs = resabs+w21b(k)*(dabs(fval1)+dabs(fval2))
savfun(ipx) = fval
fv3(k) = fval1
fv4(k) = fval2
15 continue
c
c test for convergence.
c
result = res21*hlgth
resabs = resabs*dhlgth
reskh = 0.5d+00*res21
resasc = w21b(6)*dabs(fcentr-reskh)
do 20 k = 1,5
resasc = resasc+w21a(k)*(dabs(fv1(k)-reskh)+dabs(fv2(k)-reskh))
* +w21b(k)*(dabs(fv3(k)-reskh)+dabs(fv4(k)-reskh))
20 continue
abserr = dabs((res21-res10)*hlgth)
resasc = resasc*dhlgth
go to 65
c
c compute the integral using the 43-point formula.
c
25 res43 = w43b(12)*fcentr
neval = 43
do 30 k=1,10
res43 = res43+savfun(k)*w43a(k)
30 continue
do 40 k=1,11
ipx = ipx+1
absc = hlgth*x3(k)
fval = f(absc+centr)+f(centr-absc)
res43 = res43+fval*w43b(k)
savfun(ipx) = fval
40 continue
c
c test for convergence.
c
result = res43*hlgth
abserr = dabs((res43-res21)*hlgth)
go to 65
c
c compute the integral using the 87-point formula.
c
45 res87 = w87b(23)*fcentr
neval = 87
do 50 k=1,21
res87 = res87+savfun(k)*w87a(k)
50 continue
do 60 k=1,22
absc = hlgth*x4(k)
res87 = res87+w87b(k)*(f(absc+centr)+f(centr-absc))
60 continue
result = res87*hlgth
abserr = dabs((res87-res43)*hlgth)
65 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)
if (abserr.le.dmax1(epsabs,epsrel*dabs(result))) ier = 0
c ***jump out of do-loop
if (ier.eq.0) go to 999
70 continue
80 call xerror('abnormal return from dqng ',26,ier,0)
999 return
end
|
