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
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