dqk51.f
9.32 KB
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
!
! fvn comment :
! Unmodified quadpack routine from http://www.netlib.org/quadpack
!
subroutine dqk51(f,a,b,result,abserr,resabs,resasc)
!***begin prologue dqk51
!***date written 800101 (yymmdd)
!***revision date 830518 (yymmdd)
!***category no. h2a1a2
!***keywords 51-point gauss-kronrod rules
!***author piessens,robert,appl. math. & progr. div. - k.u.leuven
! de doncker,elise,appl. math & progr. div. - k.u.leuven
!***purpose to compute i = integral of f over (a,b) with error
! estimate
! j = integral of abs(f) over (a,b)
!***description
!
! integration rules
! standard fortran subroutine
! double precision version
!
! parameters
! on entry
! f - double precision
! function subroutine defining the integrand
! function f(x). the actual name for f needs to be
! declared e x t e r n a l in the calling program.
!
! a - double precision
! lower limit of integration
!
! b - double precision
! upper limit of integration
!
! on return
! result - double precision
! approximation to the integral i
! result is computed by applying the 51-point
! kronrod rule (resk) obtained by optimal addition
! of abscissae to the 25-point gauss rule (resg).
!
! abserr - double precision
! estimate of the modulus of the absolute error,
! which should not exceed abs(i-result)
!
! resabs - double precision
! approximation to the integral j
!
! resasc - double precision
! approximation to the integral of abs(f-i/(b-a))
! over (a,b)
!
!***references (none)
!***routines called d1mach
!***end prologue dqk51
!
double precision a,absc,abserr,b,centr,dabs,dhlgth,dmax1,dmin1, &
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
!
dimension fv1(25),fv2(25),xgk(26),wgk(26),wg(13)
!
! the abscissae and weights are given for the interval (-1,1).
! because of symmetry only the positive abscissae and their
! corresponding weights are given.
!
! xgk - abscissae of the 51-point kronrod rule
! xgk(2), xgk(4), ... abscissae of the 25-point
! gauss rule
! xgk(1), xgk(3), ... abscissae which are optimally
! added to the 25-point gauss rule
!
! wgk - weights of the 51-point kronrod rule
!
! wg - weights of the 25-point gauss rule
!
!
! gauss quadrature weights and kronron quadrature abscissae and weights
! as evaluated with 80 decimal digit arithmetic by l. w. fullerton,
! bell labs, nov. 1981.
!
data wg ( 1) / 0.011393798501026287947902964113235d0 /
data wg ( 2) / 0.026354986615032137261901815295299d0 /
data wg ( 3) / 0.040939156701306312655623487711646d0 /
data wg ( 4) / 0.054904695975835191925936891540473d0 /
data wg ( 5) / 0.068038333812356917207187185656708d0 /
data wg ( 6) / 0.080140700335001018013234959669111d0 /
data wg ( 7) / 0.091028261982963649811497220702892d0 /
data wg ( 8) / 0.100535949067050644202206890392686d0 /
data wg ( 9) / 0.108519624474263653116093957050117d0 /
data wg ( 10) / 0.114858259145711648339325545869556d0 /
data wg ( 11) / 0.119455763535784772228178126512901d0 /
data wg ( 12) / 0.122242442990310041688959518945852d0 /
data wg ( 13) / 0.123176053726715451203902873079050d0 /
!
data xgk ( 1) / 0.999262104992609834193457486540341d0 /
data xgk ( 2) / 0.995556969790498097908784946893902d0 /
data xgk ( 3) / 0.988035794534077247637331014577406d0 /
data xgk ( 4) / 0.976663921459517511498315386479594d0 /
data xgk ( 5) / 0.961614986425842512418130033660167d0 /
data xgk ( 6) / 0.942974571228974339414011169658471d0 /
data xgk ( 7) / 0.920747115281701561746346084546331d0 /
data xgk ( 8) / 0.894991997878275368851042006782805d0 /
data xgk ( 9) / 0.865847065293275595448996969588340d0 /
data xgk ( 10) / 0.833442628760834001421021108693570d0 /
data xgk ( 11) / 0.797873797998500059410410904994307d0 /
data xgk ( 12) / 0.759259263037357630577282865204361d0 /
data xgk ( 13) / 0.717766406813084388186654079773298d0 /
data xgk ( 14) / 0.673566368473468364485120633247622d0 /
data xgk ( 15) / 0.626810099010317412788122681624518d0 /
data xgk ( 16) / 0.577662930241222967723689841612654d0 /
data xgk ( 17) / 0.526325284334719182599623778158010d0 /
data xgk ( 18) / 0.473002731445714960522182115009192d0 /
data xgk ( 19) / 0.417885382193037748851814394594572d0 /
data xgk ( 20) / 0.361172305809387837735821730127641d0 /
data xgk ( 21) / 0.303089538931107830167478909980339d0 /
data xgk ( 22) / 0.243866883720988432045190362797452d0 /
data xgk ( 23) / 0.183718939421048892015969888759528d0 /
data xgk ( 24) / 0.122864692610710396387359818808037d0 /
data xgk ( 25) / 0.061544483005685078886546392366797d0 /
data xgk ( 26) / 0.000000000000000000000000000000000d0 /
!
data wgk ( 1) / 0.001987383892330315926507851882843d0 /
data wgk ( 2) / 0.005561932135356713758040236901066d0 /
data wgk ( 3) / 0.009473973386174151607207710523655d0 /
data wgk ( 4) / 0.013236229195571674813656405846976d0 /
data wgk ( 5) / 0.016847817709128298231516667536336d0 /
data wgk ( 6) / 0.020435371145882835456568292235939d0 /
data wgk ( 7) / 0.024009945606953216220092489164881d0 /
data wgk ( 8) / 0.027475317587851737802948455517811d0 /
data wgk ( 9) / 0.030792300167387488891109020215229d0 /
data wgk ( 10) / 0.034002130274329337836748795229551d0 /
data wgk ( 11) / 0.037116271483415543560330625367620d0 /
data wgk ( 12) / 0.040083825504032382074839284467076d0 /
data wgk ( 13) / 0.042872845020170049476895792439495d0 /
data wgk ( 14) / 0.045502913049921788909870584752660d0 /
data wgk ( 15) / 0.047982537138836713906392255756915d0 /
data wgk ( 16) / 0.050277679080715671963325259433440d0 /
data wgk ( 17) / 0.052362885806407475864366712137873d0 /
data wgk ( 18) / 0.054251129888545490144543370459876d0 /
data wgk ( 19) / 0.055950811220412317308240686382747d0 /
data wgk ( 20) / 0.057437116361567832853582693939506d0 /
data wgk ( 21) / 0.058689680022394207961974175856788d0 /
data wgk ( 22) / 0.059720340324174059979099291932562d0 /
data wgk ( 23) / 0.060539455376045862945360267517565d0 /
data wgk ( 24) / 0.061128509717053048305859030416293d0 /
data wgk ( 25) / 0.061471189871425316661544131965264d0 /
! note: wgk (26) was calculated from the values of wgk(1..25)
data wgk ( 26) / 0.061580818067832935078759824240066d0 /
!
!
! list of major variables
! -----------------------
!
! centr - mid point of the interval
! hlgth - half-length of the interval
! absc - abscissa
! fval* - function value
! resg - result of the 25-point gauss formula
! resk - result of the 51-point kronrod formula
! reskh - approximation to the mean value of f over (a,b),
! i.e. to i/(b-a)
!
! machine dependent constants
! ---------------------------
!
! epmach is the largest relative spacing.
! uflow is the smallest positive magnitude.
!
!***first executable statement dqk51
epmach = d1mach(4)
uflow = d1mach(1)
!
centr = 0.5d+00*(a+b)
hlgth = 0.5d+00*(b-a)
dhlgth = dabs(hlgth)
!
! compute the 51-point kronrod approximation to
! the integral, and estimate the absolute error.
!
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 subroutine