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
|
subroutine dqk51(f,a,b,result,abserr,resabs,resasc)
c***begin prologue dqk51
c***date written 800101 (yymmdd)
c***revision date 830518 (yymmdd)
c***category no. h2a1a2
c***keywords 51-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 subroutine 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 51-point
c kronrod rule (resk) obtained by optimal addition
c of abscissae to the 25-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 dqk51
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(25),fv2(25),xgk(26),wgk(26),wg(13)
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 51-point kronrod rule
c xgk(2), xgk(4), ... abscissae of the 25-point
c gauss rule
c xgk(1), xgk(3), ... abscissae which are optimally
c added to the 25-point gauss rule
c
c wgk - weights of the 51-point kronrod rule
c
c wg - weights of the 25-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.0113937985 0102628794 7902964113 235 d0 /
data wg ( 2) / 0.0263549866 1503213726 1901815295 299 d0 /
data wg ( 3) / 0.0409391567 0130631265 5623487711 646 d0 /
data wg ( 4) / 0.0549046959 7583519192 5936891540 473 d0 /
data wg ( 5) / 0.0680383338 1235691720 7187185656 708 d0 /
data wg ( 6) / 0.0801407003 3500101801 3234959669 111 d0 /
data wg ( 7) / 0.0910282619 8296364981 1497220702 892 d0 /
data wg ( 8) / 0.1005359490 6705064420 2206890392 686 d0 /
data wg ( 9) / 0.1085196244 7426365311 6093957050 117 d0 /
data wg ( 10) / 0.1148582591 4571164833 9325545869 556 d0 /
data wg ( 11) / 0.1194557635 3578477222 8178126512 901 d0 /
data wg ( 12) / 0.1222424429 9031004168 8959518945 852 d0 /
data wg ( 13) / 0.1231760537 2671545120 3902873079 050 d0 /
c
data xgk ( 1) / 0.9992621049 9260983419 3457486540 341 d0 /
data xgk ( 2) / 0.9955569697 9049809790 8784946893 902 d0 /
data xgk ( 3) / 0.9880357945 3407724763 7331014577 406 d0 /
data xgk ( 4) / 0.9766639214 5951751149 8315386479 594 d0 /
data xgk ( 5) / 0.9616149864 2584251241 8130033660 167 d0 /
data xgk ( 6) / 0.9429745712 2897433941 4011169658 471 d0 /
data xgk ( 7) / 0.9207471152 8170156174 6346084546 331 d0 /
data xgk ( 8) / 0.8949919978 7827536885 1042006782 805 d0 /
data xgk ( 9) / 0.8658470652 9327559544 8996969588 340 d0 /
data xgk ( 10) / 0.8334426287 6083400142 1021108693 570 d0 /
data xgk ( 11) / 0.7978737979 9850005941 0410904994 307 d0 /
data xgk ( 12) / 0.7592592630 3735763057 7282865204 361 d0 /
data xgk ( 13) / 0.7177664068 1308438818 6654079773 298 d0 /
data xgk ( 14) / 0.6735663684 7346836448 5120633247 622 d0 /
data xgk ( 15) / 0.6268100990 1031741278 8122681624 518 d0 /
data xgk ( 16) / 0.5776629302 4122296772 3689841612 654 d0 /
data xgk ( 17) / 0.5263252843 3471918259 9623778158 010 d0 /
data xgk ( 18) / 0.4730027314 4571496052 2182115009 192 d0 /
data xgk ( 19) / 0.4178853821 9303774885 1814394594 572 d0 /
data xgk ( 20) / 0.3611723058 0938783773 5821730127 641 d0 /
data xgk ( 21) / 0.3030895389 3110783016 7478909980 339 d0 /
data xgk ( 22) / 0.2438668837 2098843204 5190362797 452 d0 /
data xgk ( 23) / 0.1837189394 2104889201 5969888759 528 d0 /
data xgk ( 24) / 0.1228646926 1071039638 7359818808 037 d0 /
data xgk ( 25) / 0.0615444830 0568507888 6546392366 797 d0 /
data xgk ( 26) / 0.0000000000 0000000000 0000000000 000 d0 /
c
data wgk ( 1) / 0.0019873838 9233031592 6507851882 843 d0 /
data wgk ( 2) / 0.0055619321 3535671375 8040236901 066 d0 /
data wgk ( 3) / 0.0094739733 8617415160 7207710523 655 d0 /
data wgk ( 4) / 0.0132362291 9557167481 3656405846 976 d0 /
data wgk ( 5) / 0.0168478177 0912829823 1516667536 336 d0 /
data wgk ( 6) / 0.0204353711 4588283545 6568292235 939 d0 /
data wgk ( 7) / 0.0240099456 0695321622 0092489164 881 d0 /
data wgk ( 8) / 0.0274753175 8785173780 2948455517 811 d0 /
data wgk ( 9) / 0.0307923001 6738748889 1109020215 229 d0 /
data wgk ( 10) / 0.0340021302 7432933783 6748795229 551 d0 /
data wgk ( 11) / 0.0371162714 8341554356 0330625367 620 d0 /
data wgk ( 12) / 0.0400838255 0403238207 4839284467 076 d0 /
data wgk ( 13) / 0.0428728450 2017004947 6895792439 495 d0 /
data wgk ( 14) / 0.0455029130 4992178890 9870584752 660 d0 /
data wgk ( 15) / 0.0479825371 3883671390 6392255756 915 d0 /
data wgk ( 16) / 0.0502776790 8071567196 3325259433 440 d0 /
data wgk ( 17) / 0.0523628858 0640747586 4366712137 873 d0 /
data wgk ( 18) / 0.0542511298 8854549014 4543370459 876 d0 /
data wgk ( 19) / 0.0559508112 2041231730 8240686382 747 d0 /
data wgk ( 20) / 0.0574371163 6156783285 3582693939 506 d0 /
data wgk ( 21) / 0.0586896800 2239420796 1974175856 788 d0 /
data wgk ( 22) / 0.0597203403 2417405997 9099291932 562 d0 /
data wgk ( 23) / 0.0605394553 7604586294 5360267517 565 d0 /
data wgk ( 24) / 0.0611285097 1705304830 5859030416 293 d0 /
data wgk ( 25) / 0.0614711898 7142531666 1544131965 264 d0 /
c note: wgk (26) was calculated from the values of wgk(1..25)
data wgk ( 26) / 0.0615808180 6783293507 8759824240 066 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 25-point gauss formula
c resk - result of the 51-point kronrod formula
c reskh - approximation to the mean value of f over (a,b),
c i.e. 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 dqk51
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 51-point kronrod approximation to
c the integral, and estimate the absolute error.
c
fc = f(centr)
resg = wg(13)*fc
resk = wgk(26)*fc
resabs = dabs(resk)
do 10 j=1,12
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,13
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(26)*dabs(fc-reskh)
do 20 j=1,25
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
|
