|
2329
|
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 |
|
3135
|
20 C F - SUBROUTINE F(X,IERR,RESULT) DEFINING THE INTEGRAND |
|
2329
|
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, |
|
3136
|
53 * D1MACH,EPMACH,FC,FSUM,FVAL1,FVAL2,FV1,FV2,HLGTH,RESABS,RESASC, |
|
2329
|
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 |
|
3136
|
142 CALL F(CENTR,IERR,FC) |
|
2329
|
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) |
|
3135
|
149 CALL F(CENTR-ABSC,IERR,FVAL1) |
|
2329
|
150 IF (IERR .LT. 0) RETURN |
|
3135
|
151 CALL F(CENTR+ABSC,IERR,FVAL2) |
|
2329
|
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) |
|
3135
|
163 CALL F(CENTR-ABSC,IERR,FVAL1) |
|
2329
|
164 IF (IERR .LT. 0) RETURN |
|
3135
|
165 CALL F(CENTR+ABSC,IERR,FVAL2) |
|
2329
|
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 |
