GNU Octave  4.2.1 A high-level interpreted language, primarily intended for numerical computations, mostly compatible with Matlab
qagi.f
Go to the documentation of this file.
1  subroutine qagi(f,bound,inf,epsabs,epsrel,result,abserr,neval,
2  * ier,limit,lenw,last,iwork,work)
3 c***begin prologue qagi
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 parameters
24 c on entry
25 c f - subroutine f(x,result) defining the integrand
26 c function f(x). the actual name for f needs to be
27 c declared e x t e r n a l in the driver program.
28 c
29 c bound - real
30 c finite bound of integration range
31 c (has no meaning if interval is doubly-infinite)
32 c
33 c inf - integer
34 c indicating the kind of integration range involved
35 c inf = 1 corresponds to (bound,+infinity),
36 c inf = -1 to (-infinity,bound),
37 c inf = 2 to (-infinity,+infinity).
38 c
39 c epsabs - real
40 c absolute accuracy requested
41 c epsrel - real
42 c relative accuracy requested
43 c if epsabs.le.0
44 c and epsrel.lt.max(50*rel.mach.acc.,0.5d-28),
45 c the routine will end with ier = 6.
46 c
47 c
48 c on return
49 c result - real
50 c approximation to the integral
51 c
52 c abserr - real
53 c estimate of the modulus of the absolute error,
54 c which should equal or exceed abs(i-result)
55 c
56 c neval - integer
57 c number of integrand evaluations
58 c
59 c ier - integer
60 c ier = 0 normal and reliable termination of the
61 c routine. it is assumed that the requested
62 c accuracy has been achieved.
63 c - ier.gt.0 abnormal termination of the routine. the
64 c estimates for result and error are less
65 c reliable. it is assumed that the requested
66 c accuracy has not been achieved.
67 c error messages
68 c ier = 1 maximum number of subdivisions allowed
69 c has been achieved. one can allow more
70 c subdivisions by increasing the value of
71 c limit (and taking the according dimension
72 c adjustments into account). however, if
73 c this yields no improvement it is advised
74 c to analyze the integrand in order to
75 c determine the integration difficulties. if
76 c the position of a local difficulty can be
77 c determined (e.g. singularity,
78 c discontinuity within the interval) one
79 c will probably gain from splitting up the
80 c interval at this point and calling the
81 c integrator on the subranges. if possible,
82 c an appropriate special-purpose integrator
83 c should be used, which is designed for
84 c handling the type of difficulty involved.
85 c = 2 the occurrence of roundoff error is
86 c detected, which prevents the requested
87 c tolerance from being achieved.
88 c the error may be under-estimated.
89 c = 3 extremely bad integrand behaviour occurs
90 c at some points of the integration
91 c interval.
92 c = 4 the algorithm does not converge.
93 c roundoff error is detected in the
94 c extrapolation table.
95 c it is assumed that the requested tolerance
96 c cannot be achieved, and that the returned
97 c result is the best which can be obtained.
98 c = 5 the integral is probably divergent, or
99 c slowly convergent. it must be noted that
100 c divergence can occur with any other value
101 c of ier.
102 c = 6 the input is invalid, because
103 c (epsabs.le.0 and
104 c epsrel.lt.max(50*rel.mach.acc.,0.5d-28))
105 c or limit.lt.1 or leniw.lt.limit*4.
106 c result, abserr, neval, last are set to
107 c zero. exept when limit or leniw is
108 c invalid, iwork(1), work(limit*2+1) and
109 c work(limit*3+1) are set to zero, work(1)
110 c is set to a and work(limit+1) to b.
111 c
112 c dimensioning parameters
113 c limit - integer
114 c dimensioning parameter for iwork
115 c limit determines the maximum number of subintervals
116 c in the partition of the given integration interval
117 c (a,b), limit.ge.1.
118 c if limit.lt.1, the routine will end with ier = 6.
119 c
120 c lenw - integer
121 c dimensioning parameter for work
122 c lenw must be at least limit*4.
123 c if lenw.lt.limit*4, the routine will end
124 c with ier = 6.
125 c
126 c last - integer
127 c on return, last equals the number of subintervals
128 c produced in the subdivision process, which
129 c determines the number of significant elements
130 c actually in the work arrays.
131 c
132 c work arrays
133 c iwork - integer
134 c vector of dimension at least limit, the first
135 c k elements of which contain pointers
136 c to the error estimates over the subintervals,
137 c such that work(limit*3+iwork(1)),... ,
138 c work(limit*3+iwork(k)) form a decreasing
139 c sequence, with k = last if last.le.(limit/2+2), and
140 c k = limit+1-last otherwise
141 c
142 c work - real
143 c vector of dimension at least lenw
144 c on return
145 c work(1), ..., work(last) contain the left
146 c end points of the subintervals in the
147 c partition of (a,b),
148 c work(limit+1), ..., work(limit+last) contain
149 c the right end points,
150 c work(limit*2+1), ...,work(limit*2+last) contain the
151 c integral approximations over the subintervals,
152 c work(limit*3+1), ..., work(limit*3)
153 c contain the error estimates.
154 c***references (none)
155 c***routines called qagie,xerror
156 c***end prologue qagi
157 c
158  real abserr, epsabs,epsrel,result,work
159  integer ier,iwork, lenw,limit,lvl,l1,l2,l3,neval
160 c
161  dimension iwork(limit),work(lenw)
162 c
163  external f
164 c
165 c check validity of limit and lenw.
166 c
167 c***first executable statement qagi
168  ier = 6
169  neval = 0
170  last = 0
171  result = 0.0e+00
172  abserr = 0.0e+00
173  if(limit.lt.1.or.lenw.lt.limit*4) go to 10
174 c
175 c prepare call for qagie.
176 c
177  l1 = limit+1
178  l2 = limit+l1
179  l3 = limit+l2
180 c
181  call qagie(f,bound,inf,epsabs,epsrel,limit,result,abserr,
182  * neval,ier,work(1),work(l1),work(l2),work(l3),iwork,last)
183 c
184 c call error handler if necessary.
185 c
186  lvl = 0
187 10 if(ier.eq.6) lvl = 1
188  if(ier.ne.0) call xerror('abnormal return from qagi',26,ier,lvl)
189  return
190  end
subroutine xerror(MESSG, NMESSG, NERR, LEVEL)
Definition: xerror.f:1
subroutine qagi(f, bound, inf, epsabs, epsrel, result, abserr, neval, ier, limit, lenw, last, iwork, work)
Definition: qagi.f:1
may be zero for pure relative error test tem the relative tolerance must be greater than or equal to