GNU Octave  9.1.0
A high-level interpreted language, primarily intended for numerical computations, mostly compatible with Matlab
qk15i.f
Go to the documentation of this file.
1  subroutine qk15i(f,boun,inf,a,b,result,abserr,resabs,resasc,ierr)
2 c***begin prologue qk15i
3 c***date written 800101 (yymmdd)
4 c***revision date 830518 (yymmdd)
5 c***category no. h2a3a2,h2a4a2
6 c***keywords 15-point transformed 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 the original (infinite integration range is mapped
10 c onto the interval (0,1) and (a,b) is a part of (0,1).
11 c it is the purpose to compute
12 c i = integral of transformed integrand over (a,b),
13 c j = integral of abs(transformed integrand) over (a,b).
14 c***description
15 c
16 c integration rule
17 c standard fortran subroutine
18 c real version
19 c
20 c parameters
21 c on entry
22 c f - subroutine f(x,ierr,result) defining the integrand
23 c function f(x). the actual name for f needs to be
24 c declared e x t e r n a l in the calling program.
25 c
26 c boun - real
27 c finite bound of original integration
28 c range (set to zero if inf = +2)
29 c
30 c inf - integer
31 c if inf = -1, the original interval is
32 c (-infinity,bound),
33 c if inf = +1, the original interval is
34 c (bound,+infinity),
35 c if inf = +2, the original interval is
36 c (-infinity,+infinity) and
37 c the integral is computed as the sum of two
38 c integrals, one over (-infinity,0) and one over
39 c (0,+infinity).
40 c
41 c a - real
42 c lower limit for integration over subrange
43 c of (0,1)
44 c
45 c b - real
46 c upper limit for integration over subrange
47 c of (0,1)
48 c
49 c on return
50 c result - real
51 c approximation to the integral i
52 c result is computed by applying the 15-point
53 c kronrod rule(resk) obtained by optimal addition
54 c of abscissae to the 7-point gauss rule(resg).
55 c
56 c abserr - real
57 c estimate of the modulus of the absolute error,
58 c which should equal or exceed abs(i-result)
59 c
60 c resabs - real
61 c approximation to the integral j
62 c
63 c resasc - real
64 c approximation to the integral of
65 c abs((transformed integrand)-i/(b-a)) over (a,b)
66 c
67 c***references (none)
68 c***routines called r1mach
69 c***end prologue qk15i
70 c
71  real a,absc,absc1,absc2,abserr,b,boun,centr,
72  * dinf,r1mach,epmach,fc,fsum,fval1,fval2,fvalt,fv1,
73  * fv2,hlgth,resabs,resasc,resg,resk,reskh,result,tabsc1,tabsc2,
74  * uflow,wg,wgk,xgk
75  integer inf,j,min0
76  external f
77 c
78  dimension fv1(7),fv2(7),xgk(8),wgk(8),wg(8)
79 c
80 c the abscissae and weights are supplied for the interval
81 c (-1,1). because of symmetry only the positive abscissae and
82 c their corresponding weights are given.
83 c
84 c xgk - abscissae of the 15-point kronrod rule
85 c xgk(2), xgk(4), ... abscissae of the 7-point
86 c gauss rule
87 c xgk(1), xgk(3), ... abscissae which are optimally
88 c added to the 7-point gauss rule
89 c
90 c wgk - weights of the 15-point kronrod rule
91 c
92 c wg - weights of the 7-point gauss rule, corresponding
93 c to the abscissae xgk(2), xgk(4), ...
94 c wg(1), wg(3), ... are set to zero.
95 c
96  data xgk(1),xgk(2),xgk(3),xgk(4),xgk(5),xgk(6),xgk(7),
97  * xgk(8)/
98  * 0.9914553711208126e+00, 0.9491079123427585e+00,
99  * 0.8648644233597691e+00, 0.7415311855993944e+00,
100  * 0.5860872354676911e+00, 0.4058451513773972e+00,
101  * 0.2077849550078985e+00, 0.0000000000000000e+00/
102 c
103  data wgk(1),wgk(2),wgk(3),wgk(4),wgk(5),wgk(6),wgk(7),
104  * wgk(8)/
105  * 0.2293532201052922e-01, 0.6309209262997855e-01,
106  * 0.1047900103222502e+00, 0.1406532597155259e+00,
107  * 0.1690047266392679e+00, 0.1903505780647854e+00,
108  * 0.2044329400752989e+00, 0.2094821410847278e+00/
109 c
110  data wg(1),wg(2),wg(3),wg(4),wg(5),wg(6),wg(7),wg(8)/
111  * 0.0000000000000000e+00, 0.1294849661688697e+00,
112  * 0.0000000000000000e+00, 0.2797053914892767e+00,
113  * 0.0000000000000000e+00, 0.3818300505051189e+00,
114  * 0.0000000000000000e+00, 0.4179591836734694e+00/
115 c
116 c
117 c list of major variables
118 c -----------------------
119 c
120 c centr - mid point of the interval
121 c hlgth - half-length of the interval
122 c absc* - abscissa
123 c tabsc* - transformed abscissa
124 c fval* - function value
125 c resg - result of the 7-point gauss formula
126 c resk - result of the 15-point kronrod formula
127 c reskh - approximation to the mean value of the transformed
128 c integrand over (a,b), i.e. to i/(b-a)
129 c
130 c machine dependent constants
131 c ---------------------------
132 c
133 c epmach is the largest relative spacing.
134 c uflow is the smallest positive magnitude.
135 c
136 c***first executable statement qk15i
137  epmach = r1mach(4)
138  uflow = r1mach(1)
139  dinf = min0(1,inf)
140 c
141  centr = 0.5e+00*(a+b)
142  hlgth = 0.5e+00*(b-a)
143  tabsc1 = boun+dinf*(0.1e+01-centr)/centr
144  call f(tabsc1, ierr, fval1)
145  if (ierr.lt.0) return
146  if(inf.eq.2) then
147  call f(-tabsc1, ierr, fval1)
148  if (ierr.lt.0) return
149  fval1 = fval1 + fvalt
150  endif
151  fc = (fval1/centr)/centr
152 c
153 c compute the 15-point kronrod approximation to
154 c the integral, and estimate the error.
155 c
156  resg = wg(8)*fc
157  resk = wgk(8)*fc
158  resabs = abs(resk)
159  do 10 j=1,7
160  absc = hlgth*xgk(j)
161  absc1 = centr-absc
162  absc2 = centr+absc
163  tabsc1 = boun+dinf*(0.1e+01-absc1)/absc1
164  tabsc2 = boun+dinf*(0.1e+01-absc2)/absc2
165  call f(tabsc1, ierr, fval1)
166  if (ierr.lt.0) return
167  call f(tabsc2, ierr, fval2)
168  if (ierr.lt.0) return
169  if(inf.eq.2) then
170  call f(-tabsc1,ierr,fvalt)
171  if (ierr.lt.0) return
172  fval1 = fval1 + fvalt
173  endif
174  if(inf.eq.2) then
175  call f(-tabsc2,ierr,fvalt)
176  if (ierr.lt.0) return
177  fval2 = fval2 + fvalt
178  endif
179  fval1 = (fval1/absc1)/absc1
180  fval2 = (fval2/absc2)/absc2
181  fv1(j) = fval1
182  fv2(j) = fval2
183  fsum = fval1+fval2
184  resg = resg+wg(j)*fsum
185  resk = resk+wgk(j)*fsum
186  resabs = resabs+wgk(j)*(abs(fval1)+abs(fval2))
187  10 continue
188  reskh = resk*0.5e+00
189  resasc = wgk(8)*abs(fc-reskh)
190  do 20 j=1,7
191  resasc = resasc+wgk(j)*(abs(fv1(j)-reskh)+abs(fv2(j)-reskh))
192  20 continue
193  result = resk*hlgth
194  resasc = resasc*hlgth
195  resabs = resabs*hlgth
196  abserr = abs((resk-resg)*hlgth)
197  if(resasc.ne.0.0e+00.and.abserr.ne.0.e0) abserr = resasc*
198  * amin1(0.1e+01,(0.2e+03*abserr/resasc)**1.5e+00)
199  if(resabs.gt.uflow/(0.5e+02*epmach)) abserr = amax1
200  * ((epmach*0.5e+02)*resabs,abserr)
201  return
202  end
subroutine qk15i(f, boun, inf, a, b, result, abserr, resabs, resasc, ierr)
Definition: qk15i.f:2