GNU Octave  4.4.1
A high-level interpreted language, primarily intended for numerical computations, mostly compatible with Matlab
dqk21.f
Go to the documentation of this file.
1  SUBROUTINE dqk21(F,A,B,RESULT,ABSERR,RESABS,RESASC,IERR)
2 C***BEGIN PROLOGUE DQK21
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 DOUBLE PRECISION 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 - DOUBLE PRECISION
25 C LOWER LIMIT OF INTEGRATION
26 C
27 C B - DOUBLE PRECISION
28 C UPPER LIMIT OF INTEGRATION
29 C
30 C ON RETURN
31 C RESULT - DOUBLE PRECISION
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 - DOUBLE PRECISION
38 C ESTIMATE OF THE MODULUS OF THE ABSOLUTE ERROR,
39 C WHICH SHOULD NOT EXCEED ABS(I-RESULT)
40 C
41 C RESABS - DOUBLE PRECISION
42 C APPROXIMATION TO THE INTEGRAL J
43 C
44 C RESASC - DOUBLE PRECISION
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 D1MACH
50 C***END PROLOGUE DQK21
51 C
52  DOUBLE PRECISION A,ABSC,ABSERR,B,CENTR,DABS,DHLGTH,DMAX1,DMIN1,
53  * D1MACH,EPMACH,FC,FSUM,FVAL1,FVAL2,FV1,FV2,HLGTH,RESABS,RESASC,
54  * RESG,RESK,RESKH,RESULT,UFLOW,WG,WGK,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 C
75 C GAUSS QUADRATURE WEIGHTS AND KRONRON QUADRATURE ABSCISSAE AND WEIGHTS
76 C AS EVALUATED WITH 80 DECIMAL DIGIT ARITHMETIC BY L. W. FULLERTON,
77 C BELL LABS, NOV. 1981.
78 C
79  DATA wg( 1) / 0.0666713443 0868813759 3568809893 332 d0 /
80  DATA wg( 2) / 0.1494513491 5058059314 5776339657 697 d0 /
81  DATA wg( 3) / 0.2190863625 1598204399 5534934228 163 d0 /
82  DATA wg( 4) / 0.2692667193 0999635509 1226921569 469 d0 /
83  DATA wg( 5) / 0.2955242247 1475287017 3892994651 338 d0 /
84 C
85  DATA xgk( 1) / 0.9956571630 2580808073 5527280689 003 d0 /
86  DATA xgk( 2) / 0.9739065285 1717172007 7964012084 452 d0 /
87  DATA xgk( 3) / 0.9301574913 5570822600 1207180059 508 d0 /
88  DATA xgk( 4) / 0.8650633666 8898451073 2096688423 493 d0 /
89  DATA xgk( 5) / 0.7808177265 8641689706 3717578345 042 d0 /
90  DATA xgk( 6) / 0.6794095682 9902440623 4327365114 874 d0 /
91  DATA xgk( 7) / 0.5627571346 6860468333 9000099272 694 d0 /
92  DATA xgk( 8) / 0.4333953941 2924719079 9265943165 784 d0 /
93  DATA xgk( 9) / 0.2943928627 0146019813 1126603103 866 d0 /
94  DATA xgk( 10) / 0.1488743389 8163121088 4826001129 720 d0 /
95  DATA xgk( 11) / 0.0000000000 0000000000 0000000000 000 d0 /
96 C
97  DATA wgk( 1) / 0.0116946388 6737187427 8064396062 192 d0 /
98  DATA wgk( 2) / 0.0325581623 0796472747 8818972459 390 d0 /
99  DATA wgk( 3) / 0.0547558965 7435199603 1381300244 580 d0 /
100  DATA wgk( 4) / 0.0750396748 1091995276 7043140916 190 d0 /
101  DATA wgk( 5) / 0.0931254545 8369760553 5065465083 366 d0 /
102  DATA wgk( 6) / 0.1093871588 0229764189 9210590325 805 d0 /
103  DATA wgk( 7) / 0.1234919762 6206585107 7958109831 074 d0 /
104  DATA wgk( 8) / 0.1347092173 1147332592 8054001771 707 d0 /
105  DATA wgk( 9) / 0.1427759385 7706008079 7094273138 717 d0 /
106  DATA wgk( 10) / 0.1477391049 0133849137 4841515972 068 d0 /
107  DATA wgk( 11) / 0.1494455540 0291690566 4936468389 821 d0 /
108 C
109 C
110 C LIST OF MAJOR VARIABLES
111 C -----------------------
112 C
113 C CENTR - MID POINT OF THE INTERVAL
114 C HLGTH - HALF-LENGTH OF THE INTERVAL
115 C ABSC - ABSCISSA
116 C FVAL* - FUNCTION VALUE
117 C RESG - RESULT OF THE 10-POINT GAUSS FORMULA
118 C RESK - RESULT OF THE 21-POINT KRONROD FORMULA
119 C RESKH - APPROXIMATION TO THE MEAN VALUE OF F OVER (A,B),
120 C I.E. TO I/(B-A)
121 C
122 C
123 C MACHINE DEPENDENT CONSTANTS
124 C ---------------------------
125 C
126 C EPMACH IS THE LARGEST RELATIVE SPACING.
127 C UFLOW IS THE SMALLEST POSITIVE MAGNITUDE.
128 C
129 C***FIRST EXECUTABLE STATEMENT DQK21
130  epmach = d1mach(4)
131  uflow = d1mach(1)
132 C
133  centr = 0.5d+00*(a+b)
134  hlgth = 0.5d+00*(b-a)
135  dhlgth = dabs(hlgth)
136 C
137 C COMPUTE THE 21-POINT KRONROD APPROXIMATION TO
138 C THE INTEGRAL, AND ESTIMATE THE ABSOLUTE ERROR.
139 C
140  resg = 0.0d+00
141  ierr = 0
142  CALL f(centr,ierr,fc)
143  IF (ierr .LT. 0) RETURN
144  resk = wgk(11)*fc
145  resabs = dabs(resk)
146  DO 10 j=1,5
147  jtw = 2*j
148  absc = hlgth*xgk(jtw)
149  CALL f(centr-absc,ierr,fval1)
150  IF (ierr .LT. 0) RETURN
151  CALL f(centr+absc,ierr,fval2)
152  IF (ierr .LT. 0) RETURN
153  fv1(jtw) = fval1
154  fv2(jtw) = fval2
155  fsum = fval1+fval2
156  resg = resg+wg(j)*fsum
157  resk = resk+wgk(jtw)*fsum
158  resabs = resabs+wgk(jtw)*(dabs(fval1)+dabs(fval2))
159  10 CONTINUE
160  DO 15 j = 1,5
161  jtwm1 = 2*j-1
162  absc = hlgth*xgk(jtwm1)
163  CALL f(centr-absc,ierr,fval1)
164  IF (ierr .LT. 0) RETURN
165  CALL f(centr+absc,ierr,fval2)
166  IF (ierr .LT. 0) RETURN
167  fv1(jtwm1) = fval1
168  fv2(jtwm1) = fval2
169  fsum = fval1+fval2
170  resk = resk+wgk(jtwm1)*fsum
171  resabs = resabs+wgk(jtwm1)*(dabs(fval1)+dabs(fval2))
172  15 CONTINUE
173  reskh = resk*0.5d+00
174  resasc = wgk(11)*dabs(fc-reskh)
175  DO 20 j=1,10
176  resasc = resasc+wgk(j)*(dabs(fv1(j)-reskh)+dabs(fv2(j)-reskh))
177  20 CONTINUE
178  result = resk*hlgth
179  resabs = resabs*dhlgth
180  resasc = resasc*dhlgth
181  abserr = dabs((resk-resg)*hlgth)
182  IF(resasc.NE.0.0d+00.AND.abserr.NE.0.0d+00)
183  * abserr = resasc*dmin1(0.1d+01,(0.2d+03*abserr/resasc)**1.5d+00)
184  IF(resabs.GT.uflow/(0.5d+02*epmach)) abserr = dmax1
185  * ((epmach*0.5d+02)*resabs,abserr)
186  RETURN
187  END
dimension
OCTAVE_EXPORT octave_value_list etc The functions then dimension(columns)