1
subroutine dqk21(f,a,b,result,abserr,resabs,resasc)
2
c***begin prologue dqk21
3
c***date written 800101 (yymmdd)
4
c***revision date 830518 (yymmdd)
5
c***category no. h2a1a2
6
c***keywords 21-point gauss-kronrod rules
7
c***author piessens,robert,appl. math. & progr. div. - k.u.leuven
8
c de doncker,elise,appl. math. & progr. div. - k.u.leuven
9
c***purpose to compute i = integral of f over (a,b), with error
11
c j = integral of abs(f) over (a,b)
15
c standard fortran subroutine
16
c double precision version
20
c f - double precision
21
c function subprogram defining the integrand
22
c function f(x). the actual name for f needs to be
23
c declared e x t e r n a l in the driver program.
25
c a - double precision
26
c lower limit of integration
28
c b - double precision
29
c upper limit of integration
32
c result - double precision
33
c approximation to the integral i
34
c result is computed by applying the 21-point
35
c kronrod rule (resk) obtained by optimal addition
36
c of abscissae to the 10-point gauss rule (resg).
38
c abserr - double precision
39
c estimate of the modulus of the absolute error,
40
c which should not exceed abs(i-result)
42
c resabs - double precision
43
c approximation to the integral j
45
c resasc - double precision
46
c approximation to the integral of abs(f-i/(b-a))
50
c***routines called d1mach
51
c***end prologue dqk21
53
double precision a,absc,abserr,b,centr,dabs,dhlgth,dmax1,dmin1,
54
* d1mach,epmach,f,fc,fsum,fval1,fval2,fv1,fv2,hlgth,resabs,resasc,
55
* resg,resk,reskh,result,uflow,wg,wgk,xgk
59
dimension fv1(10),fv2(10),wg(5),wgk(11),xgk(11)
61
c the abscissae and weights are given for the interval (-1,1).
62
c because of symmetry only the positive abscissae and their
63
c corresponding weights are given.
65
c xgk - abscissae of the 21-point kronrod rule
66
c xgk(2), xgk(4), ... abscissae of the 10-point
68
c xgk(1), xgk(3), ... abscissae which are optimally
69
c added to the 10-point gauss rule
71
c wgk - weights of the 21-point kronrod rule
73
c wg - weights of the 10-point gauss rule
76
c gauss quadrature weights and kronron quadrature abscissae and weights
77
c as evaluated with 80 decimal digit arithmetic by l. w. fullerton,
78
c bell labs, nov. 1981.
80
data wg ( 1) / 0.0666713443 0868813759 3568809893 332 d0 /
81
data wg ( 2) / 0.1494513491 5058059314 5776339657 697 d0 /
82
data wg ( 3) / 0.2190863625 1598204399 5534934228 163 d0 /
83
data wg ( 4) / 0.2692667193 0999635509 1226921569 469 d0 /
84
data wg ( 5) / 0.2955242247 1475287017 3892994651 338 d0 /
86
data xgk ( 1) / 0.9956571630 2580808073 5527280689 003 d0 /
87
data xgk ( 2) / 0.9739065285 1717172007 7964012084 452 d0 /
88
data xgk ( 3) / 0.9301574913 5570822600 1207180059 508 d0 /
89
data xgk ( 4) / 0.8650633666 8898451073 2096688423 493 d0 /
90
data xgk ( 5) / 0.7808177265 8641689706 3717578345 042 d0 /
91
data xgk ( 6) / 0.6794095682 9902440623 4327365114 874 d0 /
92
data xgk ( 7) / 0.5627571346 6860468333 9000099272 694 d0 /
93
data xgk ( 8) / 0.4333953941 2924719079 9265943165 784 d0 /
94
data xgk ( 9) / 0.2943928627 0146019813 1126603103 866 d0 /
95
data xgk ( 10) / 0.1488743389 8163121088 4826001129 720 d0 /
96
data xgk ( 11) / 0.0000000000 0000000000 0000000000 000 d0 /
98
data wgk ( 1) / 0.0116946388 6737187427 8064396062 192 d0 /
99
data wgk ( 2) / 0.0325581623 0796472747 8818972459 390 d0 /
100
data wgk ( 3) / 0.0547558965 7435199603 1381300244 580 d0 /
101
data wgk ( 4) / 0.0750396748 1091995276 7043140916 190 d0 /
102
data wgk ( 5) / 0.0931254545 8369760553 5065465083 366 d0 /
103
data wgk ( 6) / 0.1093871588 0229764189 9210590325 805 d0 /
104
data wgk ( 7) / 0.1234919762 6206585107 7958109831 074 d0 /
105
data wgk ( 8) / 0.1347092173 1147332592 8054001771 707 d0 /
106
data wgk ( 9) / 0.1427759385 7706008079 7094273138 717 d0 /
107
data wgk ( 10) / 0.1477391049 0133849137 4841515972 068 d0 /
108
data wgk ( 11) / 0.1494455540 0291690566 4936468389 821 d0 /
111
c list of major variables
112
c -----------------------
114
c centr - mid point of the interval
115
c hlgth - half-length of the interval
117
c fval* - function value
118
c resg - result of the 10-point gauss formula
119
c resk - result of the 21-point kronrod formula
120
c reskh - approximation to the mean value of f over (a,b),
124
c machine dependent constants
125
c ---------------------------
127
c epmach is the largest relative spacing.
128
c uflow is the smallest positive magnitude.
130
c***first executable statement dqk21
134
centr = 0.5d+00*(a+b)
135
hlgth = 0.5d+00*(b-a)
138
c compute the 21-point kronrod approximation to
139
c the integral, and estimate the absolute error.
147
absc = hlgth*xgk(jtw)
148
fval1 = f(centr-absc)
149
fval2 = f(centr+absc)
153
resg = resg+wg(j)*fsum
154
resk = resk+wgk(jtw)*fsum
155
resabs = resabs+wgk(jtw)*(dabs(fval1)+dabs(fval2))
159
absc = hlgth*xgk(jtwm1)
160
fval1 = f(centr-absc)
161
fval2 = f(centr+absc)
165
resk = resk+wgk(jtwm1)*fsum
166
resabs = resabs+wgk(jtwm1)*(dabs(fval1)+dabs(fval2))
169
resasc = wgk(11)*dabs(fc-reskh)
171
resasc = resasc+wgk(j)*(dabs(fv1(j)-reskh)+dabs(fv2(j)-reskh))
174
resabs = resabs*dhlgth
175
resasc = resasc*dhlgth
176
abserr = dabs((resk-resg)*hlgth)
177
if(resasc.ne.0.0d+00.and.abserr.ne.0.0d+00)
178
* abserr = resasc*dmin1(0.1d+01,(0.2d+03*abserr/resasc)**1.5d+00)
179
if(resabs.gt.uflow/(0.5d+02*epmach)) abserr = dmax1
180
* ((epmach*0.5d+02)*resabs,abserr)
added
removed
Lib/integrate/quadpack/dqk21.f
