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
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
|
subroutine dqawo(f,a,b,omega,integr,epsabs,epsrel,result,abserr,
* neval,ier,leniw,maxp1,lenw,last,iwork,work)
c***begin prologue dqawo
c***date written 800101 (yymmdd)
c***revision date 830518 (yymmdd)
c***category no. h2a2a1
c***keywords automatic integrator, special-purpose,
c integrand with oscillatory cos or sin factor,
c clenshaw-curtis method, (end point) singularities,
c extrapolation, globally adaptive
c***author piessens,robert,appl. math. & progr. div. - k.u.leuven
c de doncker,elise,appl. math. & progr. div. - k.u.leuven
c***purpose the routine calculates an approximation result to a given
c definite integral i=integral of f(x)*w(x) over (a,b)
c where w(x) = cos(omega*x)
c or w(x) = sin(omega*x),
c hopefully satisfying following claim for accuracy
c abs(i-result).le.max(epsabs,epsrel*abs(i)).
c***description
c
c computation of oscillatory integrals
c standard fortran subroutine
c double precision version
c
c parameters
c on entry
c f - double precision
c function subprogram defining the function
c f(x). the actual name for f needs to be
c declared e x t e r n a l in the driver program.
c
c a - double precision
c lower limit of integration
c
c b - double precision
c upper limit of integration
c
c omega - double precision
c parameter in the integrand weight function
c
c integr - integer
c indicates which of the weight functions is used
c integr = 1 w(x) = cos(omega*x)
c integr = 2 w(x) = sin(omega*x)
c if integr.ne.1.and.integr.ne.2, the routine will
c end with ier = 6.
c
c epsabs - double precision
c absolute accuracy requested
c epsrel - double precision
c relative accuracy requested
c if epsabs.le.0 and
c epsrel.lt.max(50*rel.mach.acc.,0.5d-28),
c the routine will end with ier = 6.
c
c on return
c result - double precision
c approximation to the integral
c
c abserr - double precision
c estimate of the modulus of the absolute error,
c which should equal or exceed abs(i-result)
c
c neval - integer
c number of integrand evaluations
c
c ier - integer
c ier = 0 normal and reliable termination of the
c routine. it is assumed that the requested
c accuracy has been achieved.
c - ier.gt.0 abnormal termination of the routine.
c the estimates for integral and error are
c less reliable. it is assumed that the
c requested accuracy has not been achieved.
c error messages
c ier = 1 maximum number of subdivisions allowed
c (= leniw/2) has been achieved. one can
c allow more subdivisions by increasing the
c value of leniw (and taking the according
c dimension adjustments into account).
c however, if this yields no improvement it
c is advised to analyze the integrand in
c order to determine the integration
c difficulties. if the position of a local
c difficulty can be determined (e.g.
c singularity, discontinuity within the
c interval) one will probably gain from
c splitting up the interval at this point
c and calling the integrator on the
c subranges. if possible, an appropriate
c special-purpose integrator should be used
c which is designed for handling the type of
c difficulty involved.
c = 2 the occurrence of roundoff error is
c detected, which prevents the requested
c tolerance from being achieved.
c the error may be under-estimated.
c = 3 extremely bad integrand behaviour occurs
c at some interior points of the
c integration interval.
c = 4 the algorithm does not converge.
c roundoff error is detected in the
c extrapolation table. it is presumed that
c the requested tolerance cannot be achieved
c due to roundoff in the extrapolation
c table, and that the returned result is
c the best which can be obtained.
c = 5 the integral is probably divergent, or
c slowly convergent. it must be noted that
c divergence can occur with any other value
c of ier.
c = 6 the input is invalid, because
c (epsabs.le.0 and
c epsrel.lt.max(50*rel.mach.acc.,0.5d-28))
c or (integr.ne.1 and integr.ne.2),
c or leniw.lt.2 or maxp1.lt.1 or
c lenw.lt.leniw*2+maxp1*25.
c result, abserr, neval, last are set to
c zero. except when leniw, maxp1 or lenw are
c invalid, work(limit*2+1), work(limit*3+1),
c iwork(1), iwork(limit+1) are set to zero,
c work(1) is set to a and work(limit+1) to
c b.
c
c dimensioning parameters
c leniw - integer
c dimensioning parameter for iwork.
c leniw/2 equals the maximum number of subintervals
c allowed in the partition of the given integration
c interval (a,b), leniw.ge.2.
c if leniw.lt.2, the routine will end with ier = 6.
c
c maxp1 - integer
c gives an upper bound on the number of chebyshev
c moments which can be stored, i.e. for the
c intervals of lengths abs(b-a)*2**(-l),
c l=0,1, ..., maxp1-2, maxp1.ge.1
c if maxp1.lt.1, the routine will end with ier = 6.
c
c lenw - integer
c dimensioning parameter for work
c lenw must be at least leniw*2+maxp1*25.
c if lenw.lt.(leniw*2+maxp1*25), the routine will
c end with ier = 6.
c
c last - integer
c on return, last equals the number of subintervals
c produced in the subdivision process, which
c determines the number of significant elements
c actually in the work arrays.
c
c work arrays
c iwork - integer
c vector of dimension at least leniw
c on return, the first k elements of which contain
c pointers to the error estimates over the
c subintervals, such that work(limit*3+iwork(1)), ..
c work(limit*3+iwork(k)) form a decreasing
c sequence, with limit = lenw/2 , and k = last
c if last.le.(limit/2+2), and k = limit+1-last
c otherwise.
c furthermore, iwork(limit+1), ..., iwork(limit+
c last) indicate the subdivision levels of the
c subintervals, such that iwork(limit+i) = l means
c that the subinterval numbered i is of length
c abs(b-a)*2**(1-l).
c
c work - double precision
c vector of dimension at least lenw
c on return
c work(1), ..., work(last) contain the left
c end points of the subintervals in the
c partition of (a,b),
c work(limit+1), ..., work(limit+last) contain
c the right end points,
c work(limit*2+1), ..., work(limit*2+last) contain
c the integral approximations over the
c subintervals,
c work(limit*3+1), ..., work(limit*3+last)
c contain the error estimates.
c work(limit*4+1), ..., work(limit*4+maxp1*25)
c provide space for storing the chebyshev moments.
c note that limit = lenw/2.
c
c***references (none)
c***routines called dqawoe,xerror
c***end prologue dqawo
c
double precision a,abserr,b,epsabs,epsrel,f,omega,result,work
integer ier,integr,iwork,last,limit,lenw,leniw,lvl,l1,l2,l3,l4,
* maxp1,momcom,neval
c
dimension iwork(leniw),work(lenw)
c
external f
c
c check validity of leniw, maxp1 and lenw.
c
c***first executable statement dqawo
ier = 6
neval = 0
last = 0
result = 0.0d+00
abserr = 0.0d+00
if(leniw.lt.2.or.maxp1.lt.1.or.lenw.lt.(leniw*2+maxp1*25))
* go to 10
c
c prepare call for dqawoe
c
limit = leniw/2
l1 = limit+1
l2 = limit+l1
l3 = limit+l2
l4 = limit+l3
call dqawoe(f,a,b,omega,integr,epsabs,epsrel,limit,1,maxp1,result,
* abserr,neval,ier,last,work(1),work(l1),work(l2),work(l3),
* iwork(1),iwork(l1),momcom,work(l4))
c
c call error handler if necessary
c
lvl = 0
10 if(ier.eq.6) lvl = 0
if(ier.ne.0) call xerror('abnormal return from dqawo',26,ier,lvl)
return
end
|
