GNU Octave  3.8.0
A high-level interpreted language, primarily intended for numerical computations, mostly compatible with Matlab
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Properties Friends Macros Pages
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,
9 c globally adaptive
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
Generated on Mon Dec 30 2013 03:04:47 for GNU Octave by   doxygen 1.8.1.2