GNU Octave  3.8.0 A high-level interpreted language, primarily intended for numerical computations, mostly compatible with Matlab
qagie.f
Go to the documentation of this file.
1  subroutine qagie(f,bound,inf,epsabs,epsrel,limit,result,abserr,
2  * neval,ier,alist,blist,rlist,elist,iord,last)
3 c***begin prologue qagie
4 c***date written 800101 (yymmdd)
5 c***revision date 830518 (yymmdd)
6 c***category no. h2a3a1,h2a4a1
7 c***keywords automatic integrator, infinite intervals,
8 c general-purpose, transformation, extrapolation,
10 c***author piessens,robert,appl. math & progr. div - k.u.leuven
11 c de doncker,elise,appl. math & progr. div - k.u.leuven
12 c***purpose the routine calculates an approximation result to a given
13 c integral i = integral of f over (bound,+infinity)
14 c or i = integral of f over (-infinity,bound)
15 c or i = integral of f over (-infinity,+infinity),
16 c hopefully satisfying following claim for accuracy
17 c abs(i-result).le.max(epsabs,epsrel*abs(i))
18 c***description
19 c
20 c integration over infinite intervals
21 c standard fortran subroutine
22 c
23 c f - subroutine f(x,ierr,result) defining the integrand
24 c function f(x). the actual name for f needs to be
25 c declared e x t e r n a l in the driver program.
26 c
27 c bound - real
28 c finite bound of integration range
29 c (has no meaning if interval is doubly-infinite)
30 c
31 c inf - real
32 c indicating the kind of integration range involved
33 c inf = 1 corresponds to (bound,+infinity),
34 c inf = -1 to (-infinity,bound),
35 c inf = 2 to (-infinity,+infinity).
36 c
37 c epsabs - real
38 c absolute accuracy requested
39 c epsrel - real
40 c relative accuracy requested
41 c if epsabs.le.0
42 c and epsrel.lt.max(50*rel.mach.acc.,0.5d-28),
43 c the routine will end with ier = 6.
44 c
45 c limit - integer
46 c gives an upper bound on the number of subintervals
47 c in the partition of (a,b), limit.ge.1
48 c
49 c on return
50 c result - real
51 c approximation to the integral
52 c
53 c abserr - real
54 c estimate of the modulus of the absolute error,
55 c which should equal or exceed abs(i-result)
56 c
57 c neval - integer
58 c number of integrand evaluations
59 c
60 c ier - integer
61 c ier = 0 normal and reliable termination of the
62 c routine. it is assumed that the requested
63 c accuracy has been achieved.
64 c - ier.gt.0 abnormal termination of the routine. the
65 c estimates for result and error are less
66 c reliable. it is assumed that the requested
67 c accuracy has not been achieved.
68 c error messages
69 c ier = 1 maximum number of subdivisions allowed
70 c has been achieved. one can allow more
71 c subdivisions by increasing the value of
72 c limit (and taking the according dimension
73 c adjustments into account). however,if
74 c this yields no improvement it is advised
75 c to analyze the integrand in order to
76 c determine the integration difficulties.
77 c if the position of a local difficulty can
78 c be determined (e.g. singularity,
79 c discontinuity within the interval) one
80 c will probably gain from splitting up the
81 c interval at this point and calling the
82 c integrator on the subranges. if possible,
83 c an appropriate special-purpose integrator
84 c should be used, which is designed for
85 c handling the type of difficulty involved.
86 c = 2 the occurrence of roundoff error is
87 c detected, which prevents the requested
88 c tolerance from being achieved.
89 c the error may be under-estimated.
90 c = 3 extremely bad integrand behaviour occurs
91 c at some points of the integration
92 c interval.
93 c = 4 the algorithm does not converge.
94 c roundoff error is detected in the
95 c extrapolation table.
96 c it is assumed that the requested tolerance
97 c cannot be achieved, and that the returned
98 c result is the best which can be obtained.
99 c = 5 the integral is probably divergent, or
100 c slowly convergent. it must be noted that
101 c divergence can occur with any other value
102 c of ier.
103 c = 6 the input is invalid, because
104 c (epsabs.le.0 and
105 c epsrel.lt.max(50*rel.mach.acc.,0.5d-28),
106 c result, abserr, neval, last, rlist(1),
107 c elist(1) and iord(1) are set to zero.
108 c alist(1) and blist(1) are set to 0
109 c and 1 respectively.
110 c
111 c alist - real
112 c vector of dimension at least limit, the first
113 c last elements of which are the left
114 c end points of the subintervals in the partition
115 c of the transformed integration range (0,1).
116 c
117 c blist - real
118 c vector of dimension at least limit, the first
119 c last elements of which are the right
120 c end points of the subintervals in the partition
121 c of the transformed integration range (0,1).
122 c
123 c rlist - real
124 c vector of dimension at least limit, the first
125 c last elements of which are the integral
126 c approximations on the subintervals
127 c
128 c elist - real
129 c vector of dimension at least limit, the first
130 c last elements of which are the moduli of the
131 c absolute error estimates on the subintervals
132 c
133 c iord - integer
134 c vector of dimension limit, the first k
135 c elements of which are pointers to the
136 c error estimates over the subintervals,
137 c such that elist(iord(1)), ..., elist(iord(k))
138 c form a decreasing sequence, with k = last
139 c if last.le.(limit/2+2), and k = limit+1-last
140 c otherwise
141 c
142 c last - integer
143 c number of subintervals actually produced
144 c in the subdivision process
145 c
146 c***references (none)
147 c***routines called qelg,qk15i,qpsrt,r1mach
148 c***end prologue qagie
149 c
150  real abseps,abserr,alist,area,area1,area12,area2,a1,
151  * a2,blist,boun,bound,b1,b2,correc,defabs,defab1,defab2,
152  * dres,r1mach,elist,epmach,epsabs,epsrel,erlarg,erlast,
153  * errbnd,errmax,error1,error2,erro12,errsum,ertest,oflow,resabs,
154  * reseps,result,res3la,rlist,rlist2,small,uflow
155  integer id,ier,ierro,inf,iord,iroff1,iroff2,iroff3,jupbnd,k,ksgn,
156  * ktmin,last,limit,maxerr,neval,nres,nrmax,numrl2
157  logical extrap,noext
158 c
159  dimension alist(limit),blist(limit),elist(limit),iord(limit),
160  * res3la(3),rlist(limit),rlist2(52)
161 c
162  external f
163 c
164 c the dimension of rlist2 is determined by the value of
165 c limexp in subroutine qelg.
166 c
167 c
168 c list of major variables
169 c -----------------------
170 c
171 c alist - list of left end points of all subintervals
172 c considered up to now
173 c blist - list of right end points of all subintervals
174 c considered up to now
175 c rlist(i) - approximation to the integral over
176 c (alist(i),blist(i))
177 c rlist2 - array of dimension at least (limexp+2),
178 c containing the part of the epsilon table
179 c wich is still needed for further computations
180 c elist(i) - error estimate applying to rlist(i)
181 c maxerr - pointer to the interval with largest error
182 c estimate
183 c errmax - elist(maxerr)
184 c erlast - error on the interval currently subdivided
185 c (before that subdivision has taken place)
186 c area - sum of the integrals over the subintervals
187 c errsum - sum of the errors over the subintervals
188 c errbnd - requested accuracy max(epsabs,epsrel*
189 c abs(result))
190 c *****1 - variable for the left subinterval
191 c *****2 - variable for the right subinterval
192 c last - index for subdivision
193 c nres - number of calls to the extrapolation routine
194 c numrl2 - number of elements currently in rlist2. if an
195 c appropriate approximation to the compounded
196 c integral has been obtained, it is put in
197 c rlist2(numrl2) after numrl2 has been increased
198 c by one.
199 c small - length of the smallest interval considered up
200 c to now, multiplied by 1.5
201 c erlarg - sum of the errors over the intervals larger
202 c than the smallest interval considered up to now
203 c extrap - logical variable denoting that the routine
204 c is attempting to perform extrapolation. i.e.
205 c before subdividing the smallest interval we
206 c try to decrease the value of erlarg.
207 c noext - logical variable denoting that extrapolation
208 c is no longer allowed (true-value)
209 c
210 c machine dependent constants
211 c ---------------------------
212 c
213 c epmach is the largest relative spacing.
214 c uflow is the smallest positive magnitude.
215 c oflow is the largest positive magnitude.
216 c
217  epmach = r1mach(4)
218 c
219 c test on validity of parameters
220 c -----------------------------
221 c
222 c***first executable statement qagie
223  ier = 0
224  neval = 0
225  last = 0
226  result = 0.0e+00
227  abserr = 0.0e+00
228  alist(1) = 0.0e+00
229  blist(1) = 0.1e+01
230  rlist(1) = 0.0e+00
231  elist(1) = 0.0e+00
232  iord(1) = 0
233  if(epsabs.le.0.0e+00.and.epsrel.lt.amax1(0.5e+02*epmach,0.5e-14))
234  * ier = 6
235  if(ier.eq.6) go to 999
236 c
237 c
238 c first approximation to the integral
239 c -----------------------------------
240 c
241 c determine the interval to be mapped onto (0,1).
242 c if inf = 2 the integral is computed as i = i1+i2, where
243 c i1 = integral of f over (-infinity,0),
244 c i2 = integral of f over (0,+infinity).
245 c
246  boun = bound
247  if(inf.eq.2) boun = 0.0e+00
248  call qk15i(f,boun,inf,0.0e+00,0.1e+01,result,abserr,
249  * defabs,resabs,ier)
250  if (ier.lt.0) return
251 c
252 c test on accuracy
253 c
254  last = 1
255  rlist(1) = result
256  elist(1) = abserr
257  iord(1) = 1
258  dres = abs(result)
259  errbnd = amax1(epsabs,epsrel*dres)
260  if(abserr.le.1.0e+02*epmach*defabs.and.abserr.gt.
261  * errbnd) ier = 2
262  if(limit.eq.1) ier = 1
263  if(ier.ne.0.or.(abserr.le.errbnd.and.abserr.ne.resabs).or.
264  * abserr.eq.0.0e+00) go to 130
265 c
266 c initialization
267 c --------------
268 c
269  uflow = r1mach(1)
270  oflow = r1mach(2)
271  rlist2(1) = result
272  errmax = abserr
273  maxerr = 1
274  area = result
275  errsum = abserr
276  abserr = oflow
277  nrmax = 1
278  nres = 0
279  ktmin = 0
280  numrl2 = 2
281  extrap = .false.
282  noext = .false.
283  ierro = 0
284  iroff1 = 0
285  iroff2 = 0
286  iroff3 = 0
287  ksgn = -1
288  if(dres.ge.(0.1e+01-0.5e+02*epmach)*defabs) ksgn = 1
289 c
290 c main do-loop
291 c ------------
292 c
293  do 90 last = 2,limit
294 c
295 c bisect the subinterval with nrmax-th largest
296 c error estimate.
297 c
298  a1 = alist(maxerr)
299  b1 = 0.5e+00*(alist(maxerr)+blist(maxerr))
300  a2 = b1
301  b2 = blist(maxerr)
302  erlast = errmax
303  call qk15i(f,boun,inf,a1,b1,area1,error1,resabs,defab1,ier)
304  if (ier.lt.0) return
305  call qk15i(f,boun,inf,a2,b2,area2,error2,resabs,defab2,ier)
306  if (ier.lt.0) return
307 c
308 c improve previous approximations to integral
309 c and error and test for accuracy.
310 c
311  area12 = area1+area2
312  erro12 = error1+error2
313  errsum = errsum+erro12-errmax
314  area = area+area12-rlist(maxerr)
315  if(defab1.eq.error1.or.defab2.eq.error2)go to 15
316  if(abs(rlist(maxerr)-area12).gt.0.1e-04*abs(area12)
317  * .or.erro12.lt.0.99e+00*errmax) go to 10
318  if(extrap) iroff2 = iroff2+1
319  if(.not.extrap) iroff1 = iroff1+1
320  10 if(last.gt.10.and.erro12.gt.errmax) iroff3 = iroff3+1
321  15 rlist(maxerr) = area1
322  rlist(last) = area2
323  errbnd = amax1(epsabs,epsrel*abs(area))
324 c
325 c test for roundoff error and eventually
326 c set error flag.
327 c
328  if(iroff1+iroff2.ge.10.or.iroff3.ge.20) ier = 2
329  if(iroff2.ge.5) ierro = 3
330 c
331 c set error flag in the case that the number of
332 c subintervals equals limit.
333 c
334  if(last.eq.limit) ier = 1
335 c
336 c set error flag in the case of bad integrand behaviour
337 c at some points of the integration range.
338 c
339  if(amax1(abs(a1),abs(b2)).le.(0.1e+01+0.1e+03*epmach)*
340  * (abs(a2)+0.1e+04*uflow)) ier = 4
341 c
342 c append the newly-created intervals to the list.
343 c
344  if(error2.gt.error1) go to 20
345  alist(last) = a2
346  blist(maxerr) = b1
347  blist(last) = b2
348  elist(maxerr) = error1
349  elist(last) = error2
350  go to 30
351  20 alist(maxerr) = a2
352  alist(last) = a1
353  blist(last) = b1
354  rlist(maxerr) = area2
355  rlist(last) = area1
356  elist(maxerr) = error2
357  elist(last) = error1
358 c
359 c call subroutine qpsrt to maintain the descending ordering
360 c in the list of error estimates and select the
361 c subinterval with nrmax-th largest error estimate (to be
362 c bisected next).
363 c
364  30 call qpsrt(limit,last,maxerr,errmax,elist,iord,nrmax)
365  if(errsum.le.errbnd) go to 115
366  if(ier.ne.0) go to 100
367  if(last.eq.2) go to 80
368  if(noext) go to 90
369  erlarg = erlarg-erlast
370  if(abs(b1-a1).gt.small) erlarg = erlarg+erro12
371  if(extrap) go to 40
372 c
373 c test whether the interval to be bisected next is the
374 c smallest interval.
375 c
376  if(abs(blist(maxerr)-alist(maxerr)).gt.small) go to 90
377  extrap = .true.
378  nrmax = 2
379  40 if(ierro.eq.3.or.erlarg.le.ertest) go to 60
380 c
381 c the smallest interval has the largest error.
382 c before bisecting decrease the sum of the errors
383 c over the larger intervals (erlarg) and perform
384 c extrapolation.
385 c
386  id = nrmax
387  jupbnd = last
388  if(last.gt.(2+limit/2)) jupbnd = limit+3-last
389  do 50 k = id,jupbnd
390  maxerr = iord(nrmax)
391  errmax = elist(maxerr)
392  if(abs(blist(maxerr)-alist(maxerr)).gt.small) go to 90
393  nrmax = nrmax+1
394  50 continue
395 c
396 c perform extrapolation.
397 c
398  60 numrl2 = numrl2+1
399  rlist2(numrl2) = area
400  call qelg(numrl2,rlist2,reseps,abseps,res3la,nres)
401  ktmin = ktmin+1
402  if(ktmin.gt.5.and.abserr.lt.0.1e-02*errsum) ier = 5
403  if(abseps.ge.abserr) go to 70
404  ktmin = 0
405  abserr = abseps
406  result = reseps
407  correc = erlarg
408  ertest = amax1(epsabs,epsrel*abs(reseps))
409  if(abserr.le.ertest) go to 100
410 c
411 c prepare bisection of the smallest interval.
412 c
413  70 if(numrl2.eq.1) noext = .true.
414  if(ier.eq.5) go to 100
415  maxerr = iord(1)
416  errmax = elist(maxerr)
417  nrmax = 1
418  extrap = .false.
419  small = small*0.5e+00
420  erlarg = errsum
421  go to 90
422  80 small = 0.375e+00
423  erlarg = errsum
424  ertest = errbnd
425  rlist2(2) = area
426  90 continue
427 c
428 c set final result and error estimate.
429 c ------------------------------------
430 c
431  100 if(abserr.eq.oflow) go to 115
432  if((ier+ierro).eq.0) go to 110
433  if(ierro.eq.3) abserr = abserr+correc
434  if(ier.eq.0) ier = 3
435  if(result.ne.0.0e+00.and.area.ne.0.0e+00)go to 105
436  if(abserr.gt.errsum)go to 115
437  if(area.eq.0.0e+00) go to 130
438  go to 110
439  105 if(abserr/abs(result).gt.errsum/abs(area))go to 115
440 c
441 c test on divergence
442 c
443  110 if(ksgn.eq.(-1).and.amax1(abs(result),abs(area)).le.
444  * defabs*0.1e-01) go to 130
445  if(0.1e-01.gt.(result/area).or.(result/area).gt.0.1e+03.
446  *or.errsum.gt.abs(area)) ier = 6
447  go to 130
448 c
449 c compute global integral sum.
450 c
451  115 result = 0.0e+00
452  do 120 k = 1,last
453  result = result+rlist(k)
454  120 continue
455  abserr = errsum
456  130 neval = 30*last-15
457  if(inf.eq.2) neval = 2*neval
458  if(ier.gt.2) ier=ier-1
459  999 return
460  end