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
|
C/MEMBR ADD NAME=QUARUL,SSI=0
subroutine quarul(f, a, b, result, abserr, resabs, resasc)
c
c based on quadpack routine quarul
c ******************************************************
c
c purpose
c to compute i = integral of f over (a,b), with error
c estimate
c j = integral of abs(f) over (a,b)
c
c calling sequence
c call quarul (f,a,b,result,abserr,resabs,resasc)
c
c parameters
c f - function subprogram defining the integrand
c function f(x). the actual name for f needs
c to be declared e x t e r n a l in the
c calling program
c
c a - lower limit of integration
c
c b - upper limit of integration
c
c result - approximation to the integral i.
c result is calculated by applying
c the 21-point gauss-kronrod rule
c (resk), obtained by optimal
c addition of abscissae to the
c 10-point gauss rule (resg).
c
c abserr - estimate of the modulus of the
c absolute error, which should not
c exceed abs(i-result)
c resabs - approximation to the integral j
c
c resasc - approximation to the integral of
c abs(f-i/(b-a)) over (a,b)
c
c ******************************************************
c .. scalar arguments ..
double precision a, abserr, b, resabs, resasc, result
c .. function arguments ..
double precision f
external f
c recuperation d'erreur
integer iero
common/ierajf/iero
c ..
c .. local scalars ..
double precision absc, centre, dhlgth, epmach, fc, fsum, fval1,
* fval2,hlgth, resg, resk, reskh, uflow
integer j
c .. local arrays ..
double precision fv1(10), fv2(10), wg(10), wgk(11), xgk(11)
c .. function references ..
double precision dlamch
c ..
c
c the abscissae and weights are given for the
c interval (-1,1) . because of symmetry only the
c positive abscissae and their corresponding
c weights are given.
c xgk - abscissae of the 21-point gauss-kronrod rule
c xgk(2), xgk(4), .... abscissae of the 10-point
c gauss rule
c xgk(1), xgk(3), .... abscissae which
c are optimally added to the 10-point
c gauss rule
c wgk - weights of the 21-point gauss-kronrod rule
c wg - weights of the 10-point gauss rule,
c corresponding to the abscissae xgk(2),
c xgk(4), ... wg(1), wg(3), ... are set
c to zero.
c
data xgk(1), xgk(2), xgk(3), xgk(4), xgk(5), xgk(6), xgk(7),xgk(8)
*, xgk(9), xgk(10), xgk(11) /0.99565716302580808073552728070d+00,
*0.97390652851717172007796401210d+00,
*0.93015749135570822600120718010d+00,
*0.86506336668898451073209668840d+00,
*0.78081772658641689706371757830d+00,
*0.67940956829902440623432736510d+00,
*0.56275713466860468333900009930d+00,
*0.43339539412924719079926594320d+00,
*0.29439286270146019813112660310d+00,
*0.14887433898163121088482600110d+00,0.0d+0/
data wgk(1), wgk(2), wgk(3), wgk(4), wgk(5), wgk(6), wgk(7),wgk(8)
*, wgk(9), wgk(10), wgk(11) /0.11694638867371874278064396060d-01,
*0.32558162307964727478818972460d-01,
*0.54755896574351996031381300240d-01,
*0.75039674810919952767043140920d-01,
*0.93125454583697605535065465080d-01,
*0.10938715880229764189921059030d+00,
*0.12349197626206585107795810980d+00,
*0.13470921731147332592805400180d+00,
*0.14277593857706008079709427310d+00,
*0.14773910490133849137484151600d+00,
*0.14944555400291690566493646840d+00/
data wg(1), wg(2), wg(3), wg(4), wg(5), wg(6), wg(7), wg(8),wg(9),
* wg(10) /0.0d+0,0.66671344308688137593568809890d-01,0.0d+0,
*0.14945134915058059314577633970d+00,0.0d+0,
*0.21908636251598204399553493420d+00,0.0d+0,
*0.26926671930999635509122692160d+00,0.0d+0,
*0.29552422471475287017389299470d+00/
epmach=dlamch('p')
uflow=dlamch('u')
c
c list of major variables
c ----------------------
c centre - mid point of the interval
c hlgth - half length of the interval
c absc - abscissa
c fval* - function value
c resg - 10-point gauss formula
c resk - 21-point gauss-kronrod formula
c reskh - approximation to mean value of f over
c (a,b), i.e. to i/(b-a)
c
centre = 0.50d+00*(a+b)
hlgth = 0.50d+00*(b-a)
dhlgth = abs(hlgth)
c
c compute the 21-point gauss-kronrod approximation to
c the integral, and estimate the absolute error
c
resg = 0.0d+00
fc = f(centre)
if(iero.ne.0) return
resk = wgk(11)*fc
resabs = abs(resk)
do 20 j=1,10
absc = hlgth*xgk(j)
fval1 = f(centre-absc)
if(iero.ne.0) return
fval2 = f(centre+absc)
if(iero.ne.0) return
fv1(j) = fval1
fv2(j) = fval2
fsum = fval1 + fval2
resg = resg + wg(j)*fsum
resk = resk + wgk(j)*fsum
resabs = resabs + wgk(j)*(abs(fval1)+abs(fval2))
20 continue
reskh = resk*0.50d+00
resasc = wgk(11)*abs(fc-reskh)
do 40 j=1,10
resasc = resasc + wgk(j)*(abs(fv1(j)-reskh)+abs(fv2(j)-reskh)
* )
40 continue
result = resk*hlgth
resabs = resabs*dhlgth
resasc = resasc*dhlgth
abserr = abs((resk-resg)*hlgth)
if (resasc.ne.0.0d+00 .and. abserr.ne.0.0d+00) abserr =resasc*
*min(0.10d+01,(0.20d+03*abserr/resasc)**1.50d+0)
if (resabs.gt.uflow/(0.50d+02*epmach)) abserr =max(epmach*resabs*
*0.50d+02,abserr)
return
end
|
