GNU Octave  6.2.0
A high-level interpreted language, primarily intended for numerical computations, mostly compatible with Matlab
qk21.f
Go to the documentation of this file.
1  subroutine qk21(f,a,b,result,abserr,resabs,resasc,ierr)
2 c***begin prologue qk21
3 c***date written 800101 (yymmdd)
4 c***revision date 830518 (yymmdd)
5 c***category no. h2a1a2
6 c***keywords 21-point 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 to compute i = integral of f over (a,b), with error
10 c estimate
11 c j = integral of abs(f) over (a,b)
12 c***description
13 c
14 c integration rules
15 c standard fortran subroutine
16 c real version
17 c
18 c parameters
19 c on entry
20 c f - subroutine f(x,ierr,result) defining the integrand
21 c function f(x). the actual name for f needs to be
22 c declared e x t e r n a l in the driver program.
23 c
24 c a - real
25 c lower limit of integration
26 c
27 c b - real
28 c upper limit of integration
29 c
30 c on return
31 c result - real
32 c approximation to the integral i
33 c result is computed by applying the 21-point
34 c kronrod rule (resk) obtained by optimal addition
35 c of abscissae to the 10-point gauss rule (resg).
36 c
37 c abserr - real
38 c estimate of the modulus of the absolute error,
39 c which should not exceed abs(i-result)
40 c
41 c resabs - real
42 c approximation to the integral j
43 c
44 c resasc - real
45 c approximation to the integral of abs(f-i/(b-a))
46 c over (a,b)
47 c
48 c***references (none)
49 c***routines called r1mach
50 c***end prologue qk21
51 c
52  real a,absc,abserr,b,centr,dhlgth,epmach,fc,fsum,fval1,fval2,
53  * fv1,fv2,hlgth,resabs,resg,resk,reskh,result,r1mach,uflow,wg,wgk,
54  * xgk
55  integer j,jtw,jtwm1
56  external f
57 c
58  dimension fv1(10),fv2(10),wg(5),wgk(11),xgk(11)
59 c
60 c the abscissae and weights are given for the interval (-1,1).
61 c because of symmetry only the positive abscissae and their
62 c corresponding weights are given.
63 c
64 c xgk - abscissae of the 21-point kronrod rule
65 c xgk(2), xgk(4), ... abscissae of the 10-point
66 c gauss rule
67 c xgk(1), xgk(3), ... abscissae which are optimally
68 c added to the 10-point gauss rule
69 c
70 c wgk - weights of the 21-point kronrod rule
71 c
72 c wg - weights of the 10-point gauss rule
73 c
74  data xgk(1),xgk(2),xgk(3),xgk(4),xgk(5),xgk(6),xgk(7),
75  * xgk(8),xgk(9),xgk(10),xgk(11)/
76  * 0.9956571630258081e+00, 0.9739065285171717e+00,
77  * 0.9301574913557082e+00, 0.8650633666889845e+00,
78  * 0.7808177265864169e+00, 0.6794095682990244e+00,
79  * 0.5627571346686047e+00, 0.4333953941292472e+00,
80  * 0.2943928627014602e+00, 0.1488743389816312e+00,
81  * 0.0000000000000000e+00/
82 c
83  data wgk(1),wgk(2),wgk(3),wgk(4),wgk(5),wgk(6),wgk(7),
84  * wgk(8),wgk(9),wgk(10),wgk(11)/
85  * 0.1169463886737187e-01, 0.3255816230796473e-01,
86  * 0.5475589657435200e-01, 0.7503967481091995e-01,
87  * 0.9312545458369761e-01, 0.1093871588022976e+00,
88  * 0.1234919762620659e+00, 0.1347092173114733e+00,
89  * 0.1427759385770601e+00, 0.1477391049013385e+00,
90  * 0.1494455540029169e+00/
91 c
92  data wg(1),wg(2),wg(3),wg(4),wg(5)/
93  * 0.6667134430868814e-01, 0.1494513491505806e+00,
94  * 0.2190863625159820e+00, 0.2692667193099964e+00,
95  * 0.2955242247147529e+00/
96 c
97 c
98 c list of major variables
99 c -----------------------
100 c
101 c centr - mid point of the interval
102 c hlgth - half-length of the interval
103 c absc - abscissa
104 c fval* - function value
105 c resg - result of the 10-point gauss formula
106 c resk - result of the 21-point kronrod formula
107 c reskh - approximation to the mean value of f over (a,b),
108 c i.e. to i/(b-a)
109 c
110 c
111 c machine dependent constants
112 c ---------------------------
113 c
114 c epmach is the largest relative spacing.
115 c uflow is the smallest positive magnitude.
116 c
117 c***first executable statement qk21
118  epmach = r1mach(4)
119  uflow = r1mach(1)
120 c
121  centr = 0.5e+00*(a+b)
122  hlgth = 0.5e+00*(b-a)
123  dhlgth = abs(hlgth)
124 c
125 c compute the 21-point kronrod approximation to
126 c the integral, and estimate the absolute error.
127 c
128  resg = 0.0e+00
129  call f(centr, ierr, fc)
130  if (ierr .lt. 0) return
131  resk = wgk(11)*fc
132  resabs = abs(resk)
133  do 10 j=1,5
134  jtw = 2*j
135  absc = hlgth*xgk(jtw)
136  call f(centr-absc,ierr,fval1)
137  if (ierr .lt. 0) return
138  call f(centr+absc,ierr,fval2)
139  if (ierr .lt. 0) return
140  fv1(jtw) = fval1
141  fv2(jtw) = fval2
142  fsum = fval1+fval2
143  resg = resg+wg(j)*fsum
144  resk = resk+wgk(jtw)*fsum
145  resabs = resabs+wgk(jtw)*(abs(fval1)+abs(fval2))
146  10 continue
147  do 15 j = 1,5
148  jtwm1 = 2*j-1
149  absc = hlgth*xgk(jtwm1)
150  call f(centr-absc,ierr,fval1)
151  if (ierr .lt. 0) return
152  call f(centr+absc,ierr,fval2)
153  if (ierr .lt. 0) return
154  fv1(jtwm1) = fval1
155  fv2(jtwm1) = fval2
156  fsum = fval1+fval2
157  resk = resk+wgk(jtwm1)*fsum
158  resabs = resabs+wgk(jtwm1)*(abs(fval1)+abs(fval2))
159  15 continue
160  reskh = resk*0.5e+00
161  resasc = wgk(11)*abs(fc-reskh)
162  do 20 j=1,10
163  resasc = resasc+wgk(j)*(abs(fv1(j)-reskh)+abs(fv2(j)-reskh))
164  20 continue
165  result = resk*hlgth
166  resabs = resabs*dhlgth
167  resasc = resasc*dhlgth
168  abserr = abs((resk-resg)*hlgth)
169  if(resasc.ne.0.0e+00.and.abserr.ne.0.0e+00)
170  * abserr = resasc*amin1(0.1e+01,
171  * (0.2e+03*abserr/resasc)**1.5e+00)
172  if(resabs.gt.uflow/(0.5e+02*epmach)) abserr = amax1
173  * ((epmach*0.5e+02)*resabs,abserr)
174  return
175  end
abs
static T abs(T x)
Definition: pr-output.cc:1678
qk21
subroutine qk21(f, a, b, result, abserr, resabs, resasc, ierr)
Definition: qk21.f:2