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
|
subroutine dqk15i(f,boun,inf,a,b,result,abserr,resabs,resasc)
c***begin prologue dqk15i
c***date written 800101 (yymmdd)
c***revision date 830518 (yymmdd)
c***category no. h2a3a2,h2a4a2
c***keywords 15-point transformed 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 the original (infinite integration range is mapped
c onto the interval (0,1) and (a,b) is a part of (0,1).
c it is the purpose to compute
c i = integral of transformed integrand over (a,b),
c j = integral of abs(transformed integrand) over (a,b).
c***description
c
c integration rule
c standard fortran subroutine
c double precision version
c
c parameters
c on entry
c f - double precision
c fuction 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 boun - double precision
c finite bound of original integration
c range (set to zero if inf = +2)
c
c inf - integer
c if inf = -1, the original interval is
c (-infinity,bound),
c if inf = +1, the original interval is
c (bound,+infinity),
c if inf = +2, the original interval is
c (-infinity,+infinity) and
c the integral is computed as the sum of two
c integrals, one over (-infinity,0) and one over
c (0,+infinity).
c
c a - double precision
c lower limit for integration over subrange
c of (0,1)
c
c b - double precision
c upper limit for integration over subrange
c of (0,1)
c
c on return
c result - double precision
c approximation to the integral i
c result is computed by applying the 15-point
c kronrod rule(resk) obtained by optimal addition
c of abscissae to the 7-point gauss rule(resg).
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 resabs - double precision
c approximation to the integral j
c
c resasc - double precision
c approximation to the integral of
c abs((transformed integrand)-i/(b-a)) over (a,b)
c
c***references (none)
c***routines called d1mach
c***end prologue dqk15i
c
double precision a,absc,absc1,absc2,abserr,b,boun,centr,dabs,dinf,
* dmax1,dmin1,d1mach,epmach,f,fc,fsum,fval1,fval2,fv1,fv2,hlgth,
* resabs,resasc,resg,resk,reskh,result,tabsc1,tabsc2,uflow,wg,wgk,
* xgk
integer inf,j
external f
c
dimension fv1(7),fv2(7),xgk(8),wgk(8),wg(8)
c
c the abscissae and weights are supplied for the interval
c (-1,1). because of symmetry only the positive abscissae and
c their corresponding weights are given.
c
c xgk - abscissae of the 15-point kronrod rule
c xgk(2), xgk(4), ... abscissae of the 7-point
c gauss rule
c xgk(1), xgk(3), ... abscissae which are optimally
c added to the 7-point gauss rule
c
c wgk - weights of the 15-point kronrod rule
c
c wg - weights of the 7-point gauss rule, corresponding
c to the abscissae xgk(2), xgk(4), ...
c wg(1), wg(3), ... are set to zero.
c
data wg(1) / 0.0d0 /
data wg(2) / 0.1294849661 6886969327 0611432679 082d0 /
data wg(3) / 0.0d0 /
data wg(4) / 0.2797053914 8927666790 1467771423 780d0 /
data wg(5) / 0.0d0 /
data wg(6) / 0.3818300505 0511894495 0369775488 975d0 /
data wg(7) / 0.0d0 /
data wg(8) / 0.4179591836 7346938775 5102040816 327d0 /
c
data xgk(1) / 0.9914553711 2081263920 6854697526 329d0 /
data xgk(2) / 0.9491079123 4275852452 6189684047 851d0 /
data xgk(3) / 0.8648644233 5976907278 9712788640 926d0 /
data xgk(4) / 0.7415311855 9939443986 3864773280 788d0 /
data xgk(5) / 0.5860872354 6769113029 4144838258 730d0 /
data xgk(6) / 0.4058451513 7739716690 6606412076 961d0 /
data xgk(7) / 0.2077849550 0789846760 0689403773 245d0 /
data xgk(8) / 0.0000000000 0000000000 0000000000 000d0 /
c
data wgk(1) / 0.0229353220 1052922496 3732008058 970d0 /
data wgk(2) / 0.0630920926 2997855329 0700663189 204d0 /
data wgk(3) / 0.1047900103 2225018383 9876322541 518d0 /
data wgk(4) / 0.1406532597 1552591874 5189590510 238d0 /
data wgk(5) / 0.1690047266 3926790282 6583426598 550d0 /
data wgk(6) / 0.1903505780 6478540991 3256402421 014d0 /
data wgk(7) / 0.2044329400 7529889241 4161999234 649d0 /
data wgk(8) / 0.2094821410 8472782801 2999174891 714d0 /
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 tabsc* - transformed abscissa
c fval* - function value
c resg - result of the 7-point gauss formula
c resk - result of the 15-point kronrod formula
c reskh - approximation to the mean value of the transformed
c integrand over (a,b), 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 dqk15i
epmach = d1mach(4)
uflow = d1mach(1)
dinf = min0(1,inf)
c
centr = 0.5d+00*(a+b)
hlgth = 0.5d+00*(b-a)
tabsc1 = boun+dinf*(0.1d+01-centr)/centr
fval1 = f(tabsc1)
if(inf.eq.2) fval1 = fval1+f(-tabsc1)
fc = (fval1/centr)/centr
c
c compute the 15-point kronrod approximation to
c the integral, and estimate the error.
c
resg = wg(8)*fc
resk = wgk(8)*fc
resabs = dabs(resk)
do 10 j=1,7
absc = hlgth*xgk(j)
absc1 = centr-absc
absc2 = centr+absc
tabsc1 = boun+dinf*(0.1d+01-absc1)/absc1
tabsc2 = boun+dinf*(0.1d+01-absc2)/absc2
fval1 = f(tabsc1)
fval2 = f(tabsc2)
if(inf.eq.2) fval1 = fval1+f(-tabsc1)
if(inf.eq.2) fval2 = fval2+f(-tabsc2)
fval1 = (fval1/absc1)/absc1
fval2 = (fval2/absc2)/absc2
fv1(j) = fval1
fv2(j) = fval2
fsum = fval1+fval2
resg = resg+wg(j)*fsum
resk = resk+wgk(j)*fsum
resabs = resabs+wgk(j)*(dabs(fval1)+dabs(fval2))
10 continue
reskh = resk*0.5d+00
resasc = wgk(8)*dabs(fc-reskh)
do 20 j=1,7
resasc = resasc+wgk(j)*(dabs(fv1(j)-reskh)+dabs(fv2(j)-reskh))
20 continue
result = resk*hlgth
resasc = resasc*hlgth
resabs = resabs*hlgth
abserr = dabs((resk-resg)*hlgth)
if(resasc.ne.0.0d+00.and.abserr.ne.0.d0) 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
|
