1*DECK DQK41 2 SUBROUTINE DQK41 (F, A, B, RESULT, ABSERR, RESABS, RESASC) 3C***BEGIN PROLOGUE DQK41 4C***PURPOSE To compute I = Integral of F over (A,B), with error 5C estimate 6C J = Integral of ABS(F) over (A,B) 7C***LIBRARY SLATEC (QUADPACK) 8C***CATEGORY H2A1A2 9C***TYPE DOUBLE PRECISION (QK41-S, DQK41-D) 10C***KEYWORDS 41-POINT GAUSS-KRONROD RULES, QUADPACK, QUADRATURE 11C***AUTHOR Piessens, Robert 12C Applied Mathematics and Programming Division 13C K. U. Leuven 14C de Doncker, Elise 15C Applied Mathematics and Programming Division 16C K. U. Leuven 17C***DESCRIPTION 18C 19C Integration rules 20C Standard fortran subroutine 21C Double precision version 22C 23C PARAMETERS 24C ON ENTRY 25C F - Double precision 26C Function subprogram defining the integrand 27C FUNCTION F(X). The actual name for F needs to be 28C declared E X T E R N A L in the calling program. 29C 30C A - Double precision 31C Lower limit of integration 32C 33C B - Double precision 34C Upper limit of integration 35C 36C ON RETURN 37C RESULT - Double precision 38C Approximation to the integral I 39C RESULT is computed by applying the 41-POINT 40C GAUSS-KRONROD RULE (RESK) obtained by optimal 41C addition of abscissae to the 20-POINT GAUSS 42C RULE (RESG). 43C 44C ABSERR - Double precision 45C Estimate of the modulus of the absolute error, 46C which should not exceed ABS(I-RESULT) 47C 48C RESABS - Double precision 49C Approximation to the integral J 50C 51C RESASC - Double precision 52C Approximation to the integral of ABS(F-I/(B-A)) 53C over (A,B) 54C 55C***REFERENCES (NONE) 56C***ROUTINES CALLED D1MACH 57C***REVISION HISTORY (YYMMDD) 58C 800101 DATE WRITTEN 59C 890531 Changed all specific intrinsics to generic. (WRB) 60C 890531 REVISION DATE from Version 3.2 61C 891214 Prologue converted to Version 4.0 format. (BAB) 62C***END PROLOGUE DQK41 63C 64 DOUBLE PRECISION A,ABSC,ABSERR,B,CENTR,DHLGTH, 65 1 D1MACH,EPMACH,F,FC,FSUM,FVAL1,FVAL2,FV1,FV2,HLGTH,RESABS,RESASC, 66 2 RESG,RESK,RESKH,RESULT,UFLOW,WG,WGK,XGK 67 INTEGER J,JTW,JTWM1 68 EXTERNAL F 69C 70 DIMENSION FV1(20),FV2(20),XGK(21),WGK(21),WG(10) 71C 72C THE ABSCISSAE AND WEIGHTS ARE GIVEN FOR THE INTERVAL (-1,1). 73C BECAUSE OF SYMMETRY ONLY THE POSITIVE ABSCISSAE AND THEIR 74C CORRESPONDING WEIGHTS ARE GIVEN. 75C 76C XGK - ABSCISSAE OF THE 41-POINT GAUSS-KRONROD RULE 77C XGK(2), XGK(4), ... ABSCISSAE OF THE 20-POINT 78C GAUSS RULE 79C XGK(1), XGK(3), ... ABSCISSAE WHICH ARE OPTIMALLY 80C ADDED TO THE 20-POINT GAUSS RULE 81C 82C WGK - WEIGHTS OF THE 41-POINT GAUSS-KRONROD RULE 83C 84C WG - WEIGHTS OF THE 20-POINT GAUSS RULE 85C 86C 87C GAUSS QUADRATURE WEIGHTS AND KRONROD QUADRATURE ABSCISSAE AND WEIGHTS 88C AS EVALUATED WITH 80 DECIMAL DIGIT ARITHMETIC BY L. W. FULLERTON, 89C BELL LABS, NOV. 1981. 90C 91 SAVE WG, XGK, WGK 92 DATA WG ( 1) / 0.0176140071 3915211831 1861962351 853 D0 / 93 DATA WG ( 2) / 0.0406014298 0038694133 1039952274 932 D0 / 94 DATA WG ( 3) / 0.0626720483 3410906356 9506535187 042 D0 / 95 DATA WG ( 4) / 0.0832767415 7670474872 4758143222 046 D0 / 96 DATA WG ( 5) / 0.1019301198 1724043503 6750135480 350 D0 / 97 DATA WG ( 6) / 0.1181945319 6151841731 2377377711 382 D0 / 98 DATA WG ( 7) / 0.1316886384 4917662689 8494499748 163 D0 / 99 DATA WG ( 8) / 0.1420961093 1838205132 9298325067 165 D0 / 100 DATA WG ( 9) / 0.1491729864 7260374678 7828737001 969 D0 / 101 DATA WG ( 10) / 0.1527533871 3072585069 8084331955 098 D0 / 102C 103 DATA XGK ( 1) / 0.9988590315 8827766383 8315576545 863 D0 / 104 DATA XGK ( 2) / 0.9931285991 8509492478 6122388471 320 D0 / 105 DATA XGK ( 3) / 0.9815078774 5025025919 3342994720 217 D0 / 106 DATA XGK ( 4) / 0.9639719272 7791379126 7666131197 277 D0 / 107 DATA XGK ( 5) / 0.9408226338 3175475351 9982722212 443 D0 / 108 DATA XGK ( 6) / 0.9122344282 5132590586 7752441203 298 D0 / 109 DATA XGK ( 7) / 0.8782768112 5228197607 7442995113 078 D0 / 110 DATA XGK ( 8) / 0.8391169718 2221882339 4529061701 521 D0 / 111 DATA XGK ( 9) / 0.7950414288 3755119835 0638833272 788 D0 / 112 DATA XGK ( 10) / 0.7463319064 6015079261 4305070355 642 D0 / 113 DATA XGK ( 11) / 0.6932376563 3475138480 5490711845 932 D0 / 114 DATA XGK ( 12) / 0.6360536807 2651502545 2836696226 286 D0 / 115 DATA XGK ( 13) / 0.5751404468 1971031534 2946036586 425 D0 / 116 DATA XGK ( 14) / 0.5108670019 5082709800 4364050955 251 D0 / 117 DATA XGK ( 15) / 0.4435931752 3872510319 9992213492 640 D0 / 118 DATA XGK ( 16) / 0.3737060887 1541956067 2548177024 927 D0 / 119 DATA XGK ( 17) / 0.3016278681 1491300432 0555356858 592 D0 / 120 DATA XGK ( 18) / 0.2277858511 4164507808 0496195368 575 D0 / 121 DATA XGK ( 19) / 0.1526054652 4092267550 5220241022 678 D0 / 122 DATA XGK ( 20) / 0.0765265211 3349733375 4640409398 838 D0 / 123 DATA XGK ( 21) / 0.0000000000 0000000000 0000000000 000 D0 / 124C 125 DATA WGK ( 1) / 0.0030735837 1852053150 1218293246 031 D0 / 126 DATA WGK ( 2) / 0.0086002698 5564294219 8661787950 102 D0 / 127 DATA WGK ( 3) / 0.0146261692 5697125298 3787960308 868 D0 / 128 DATA WGK ( 4) / 0.0203883734 6126652359 8010231432 755 D0 / 129 DATA WGK ( 5) / 0.0258821336 0495115883 4505067096 153 D0 / 130 DATA WGK ( 6) / 0.0312873067 7703279895 8543119323 801 D0 / 131 DATA WGK ( 7) / 0.0366001697 5820079803 0557240707 211 D0 / 132 DATA WGK ( 8) / 0.0416688733 2797368626 3788305936 895 D0 / 133 DATA WGK ( 9) / 0.0464348218 6749767472 0231880926 108 D0 / 134 DATA WGK ( 10) / 0.0509445739 2372869193 2707670050 345 D0 / 135 DATA WGK ( 11) / 0.0551951053 4828599474 4832372419 777 D0 / 136 DATA WGK ( 12) / 0.0591114008 8063957237 4967220648 594 D0 / 137 DATA WGK ( 13) / 0.0626532375 5478116802 5870122174 255 D0 / 138 DATA WGK ( 14) / 0.0658345971 3361842211 1563556969 398 D0 / 139 DATA WGK ( 15) / 0.0686486729 2852161934 5623411885 368 D0 / 140 DATA WGK ( 16) / 0.0710544235 5344406830 5790361723 210 D0 / 141 DATA WGK ( 17) / 0.0730306903 3278666749 5189417658 913 D0 / 142 DATA WGK ( 18) / 0.0745828754 0049918898 6581418362 488 D0 / 143 DATA WGK ( 19) / 0.0757044976 8455667465 9542775376 617 D0 / 144 DATA WGK ( 20) / 0.0763778676 7208073670 5502835038 061 D0 / 145 DATA WGK ( 21) / 0.0766007119 1799965644 5049901530 102 D0 / 146C 147C 148C LIST OF MAJOR VARIABLES 149C ----------------------- 150C 151C CENTR - MID POINT OF THE INTERVAL 152C HLGTH - HALF-LENGTH OF THE INTERVAL 153C ABSC - ABSCISSA 154C FVAL* - FUNCTION VALUE 155C RESG - RESULT OF THE 20-POINT GAUSS FORMULA 156C RESK - RESULT OF THE 41-POINT KRONROD FORMULA 157C RESKH - APPROXIMATION TO MEAN VALUE OF F OVER (A,B), I.E. 158C TO I/(B-A) 159C 160C MACHINE DEPENDENT CONSTANTS 161C --------------------------- 162C 163C EPMACH IS THE LARGEST RELATIVE SPACING. 164C UFLOW IS THE SMALLEST POSITIVE MAGNITUDE. 165C 166C***FIRST EXECUTABLE STATEMENT DQK41 167 EPMACH = D1MACH(4) 168 UFLOW = D1MACH(1) 169C 170 CENTR = 0.5D+00*(A+B) 171 HLGTH = 0.5D+00*(B-A) 172 DHLGTH = ABS(HLGTH) 173C 174C COMPUTE THE 41-POINT GAUSS-KRONROD APPROXIMATION TO 175C THE INTEGRAL, AND ESTIMATE THE ABSOLUTE ERROR. 176C 177 RESG = 0.0D+00 178 FC = F(CENTR) 179 RESK = WGK(21)*FC 180 RESABS = ABS(RESK) 181 DO 10 J=1,10 182 JTW = J*2 183 ABSC = HLGTH*XGK(JTW) 184 FVAL1 = F(CENTR-ABSC) 185 FVAL2 = F(CENTR+ABSC) 186 FV1(JTW) = FVAL1 187 FV2(JTW) = FVAL2 188 FSUM = FVAL1+FVAL2 189 RESG = RESG+WG(J)*FSUM 190 RESK = RESK+WGK(JTW)*FSUM 191 RESABS = RESABS+WGK(JTW)*(ABS(FVAL1)+ABS(FVAL2)) 192 10 CONTINUE 193 DO 15 J = 1,10 194 JTWM1 = J*2-1 195 ABSC = HLGTH*XGK(JTWM1) 196 FVAL1 = F(CENTR-ABSC) 197 FVAL2 = F(CENTR+ABSC) 198 FV1(JTWM1) = FVAL1 199 FV2(JTWM1) = FVAL2 200 FSUM = FVAL1+FVAL2 201 RESK = RESK+WGK(JTWM1)*FSUM 202 RESABS = RESABS+WGK(JTWM1)*(ABS(FVAL1)+ABS(FVAL2)) 203 15 CONTINUE 204 RESKH = RESK*0.5D+00 205 RESASC = WGK(21)*ABS(FC-RESKH) 206 DO 20 J=1,20 207 RESASC = RESASC+WGK(J)*(ABS(FV1(J)-RESKH)+ABS(FV2(J)-RESKH)) 208 20 CONTINUE 209 RESULT = RESK*HLGTH 210 RESABS = RESABS*DHLGTH 211 RESASC = RESASC*DHLGTH 212 ABSERR = ABS((RESK-RESG)*HLGTH) 213 IF(RESASC.NE.0.0D+00.AND.ABSERR.NE.0.D+00) 214 1 ABSERR = RESASC*MIN(0.1D+01,(0.2D+03*ABSERR/RESASC)**1.5D+00) 215 IF(RESABS.GT.UFLOW/(0.5D+02*EPMACH)) ABSERR = MAX 216 1 ((EPMACH*0.5D+02)*RESABS,ABSERR) 217 RETURN 218 END 219