GNU Octave  4.0.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
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
workspace_element::dimension
std::string dimension(void) const
Definition: workspace-element.h:74
qk21
subroutine qk21(f, a, b, result, abserr, resabs, resasc, ierr)
Definition: qk21.f:1
c_file_ptr_buf::uflow
int_type uflow(void)
Definition: c-file-ptr-stream.h:55
abs
T abs(T x)
Definition: pr-output.cc:3062
Generated on Thu Jun 4 2015 23:30:24 for GNU Octave by   doxygen 1.8.8