1*DECK QK41 2 SUBROUTINE QK41 (F, A, B, RESULT, ABSERR, RESABS, RESASC) 3C***BEGIN PROLOGUE QK41 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 SINGLE 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 Real version 22C 23C PARAMETERS 24C ON ENTRY 25C F - Real 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 - Real 31C Lower limit of integration 32C 33C B - Real 34C Upper limit of integration 35C 36C ON RETURN 37C RESULT - Real 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 - Real 45C Estimate of the modulus of the absolute error, 46C which should not exceed ABS(I-RESULT) 47C 48C RESABS - Real 49C Approximation to the integral J 50C 51C RESASC - Real 52C Approximation to the integral of ABS(F-I/(B-A)) 53C over (A,B) 54C 55C***REFERENCES (NONE) 56C***ROUTINES CALLED R1MACH 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 QK41 63C 64 REAL A,ABSC,ABSERR,B,CENTR,DHLGTH,EPMACH,F,FC,FSUM,FVAL1,FVAL2, 65 1 FV1,FV2,HLGTH,RESABS, 66 2 RESASC,RESG,RESK,RESKH,RESULT,R1MACH,UFLOW, 67 3 WG,WGK,XGK 68 INTEGER J,JTW,JTWM1 69 EXTERNAL F 70C 71 DIMENSION FV1(20),FV2(20),XGK(21),WGK(21),WG(10) 72C 73C THE ABSCISSAE AND WEIGHTS ARE GIVEN FOR THE INTERVAL (-1,1). 74C BECAUSE OF SYMMETRY ONLY THE POSITIVE ABSCISSAE AND THEIR 75C CORRESPONDING WEIGHTS ARE GIVEN. 76C 77C XGK - ABSCISSAE OF THE 41-POINT GAUSS-KRONROD RULE 78C XGK(2), XGK(4), ... ABSCISSAE OF THE 20-POINT 79C GAUSS RULE 80C XGK(1), XGK(3), ... ABSCISSAE WHICH ARE OPTIMALLY 81C ADDED TO THE 20-POINT GAUSS RULE 82C 83C WGK - WEIGHTS OF THE 41-POINT GAUSS-KRONROD RULE 84C 85C WG - WEIGHTS OF THE 20-POINT GAUSS RULE 86C 87 SAVE XGK, WGK, WG 88 DATA XGK(1),XGK(2),XGK(3),XGK(4),XGK(5),XGK(6),XGK(7),XGK(8), 89 1 XGK(9),XGK(10),XGK(11),XGK(12),XGK(13),XGK(14),XGK(15), 90 2 XGK(16),XGK(17),XGK(18),XGK(19),XGK(20),XGK(21)/ 91 3 0.9988590315882777E+00, 0.9931285991850949E+00, 92 4 0.9815078774502503E+00, 0.9639719272779138E+00, 93 5 0.9408226338317548E+00, 0.9122344282513259E+00, 94 6 0.8782768112522820E+00, 0.8391169718222188E+00, 95 7 0.7950414288375512E+00, 0.7463319064601508E+00, 96 8 0.6932376563347514E+00, 0.6360536807265150E+00, 97 9 0.5751404468197103E+00, 0.5108670019508271E+00, 98 1 0.4435931752387251E+00, 0.3737060887154196E+00, 99 2 0.3016278681149130E+00, 0.2277858511416451E+00, 100 3 0.1526054652409227E+00, 0.7652652113349733E-01, 101 4 0.0E+00 / 102 DATA WGK(1),WGK(2),WGK(3),WGK(4),WGK(5),WGK(6),WGK(7),WGK(8), 103 1 WGK(9),WGK(10),WGK(11),WGK(12),WGK(13),WGK(14),WGK(15),WGK(16), 104 2 WGK(17),WGK(18),WGK(19),WGK(20),WGK(21)/ 105 3 0.3073583718520532E-02, 0.8600269855642942E-02, 106 4 0.1462616925697125E-01, 0.2038837346126652E-01, 107 5 0.2588213360495116E-01, 0.3128730677703280E-01, 108 6 0.3660016975820080E-01, 0.4166887332797369E-01, 109 7 0.4643482186749767E-01, 0.5094457392372869E-01, 110 8 0.5519510534828599E-01, 0.5911140088063957E-01, 111 9 0.6265323755478117E-01, 0.6583459713361842E-01, 112 1 0.6864867292852162E-01, 0.7105442355344407E-01, 113 2 0.7303069033278667E-01, 0.7458287540049919E-01, 114 3 0.7570449768455667E-01, 0.7637786767208074E-01, 115 4 0.7660071191799966E-01/ 116 DATA WG(1),WG(2),WG(3),WG(4),WG(5),WG(6),WG(7),WG(8),WG(9),WG(10)/ 117 1 0.1761400713915212E-01, 0.4060142980038694E-01, 118 2 0.6267204833410906E-01, 0.8327674157670475E-01, 119 3 0.1019301198172404E+00, 0.1181945319615184E+00, 120 4 0.1316886384491766E+00, 0.1420961093183821E+00, 121 5 0.1491729864726037E+00, 0.1527533871307259E+00/ 122C 123C 124C LIST OF MAJOR VARIABLES 125C ----------------------- 126C 127C CENTR - MID POINT OF THE INTERVAL 128C HLGTH - HALF-LENGTH OF THE INTERVAL 129C ABSC - ABSCISSA 130C FVAL* - FUNCTION VALUE 131C RESG - RESULT OF THE 20-POINT GAUSS FORMULA 132C RESK - RESULT OF THE 41-POINT KRONROD FORMULA 133C RESKH - APPROXIMATION TO MEAN VALUE OF F OVER (A,B), I.E. 134C TO I/(B-A) 135C 136C MACHINE DEPENDENT CONSTANTS 137C --------------------------- 138C 139C EPMACH IS THE LARGEST RELATIVE SPACING. 140C UFLOW IS THE SMALLEST POSITIVE MAGNITUDE. 141C 142C***FIRST EXECUTABLE STATEMENT QK41 143 EPMACH = R1MACH(4) 144 UFLOW = R1MACH(1) 145C 146 CENTR = 0.5E+00*(A+B) 147 HLGTH = 0.5E+00*(B-A) 148 DHLGTH = ABS(HLGTH) 149C 150C COMPUTE THE 41-POINT GAUSS-KRONROD APPROXIMATION TO 151C THE INTEGRAL, AND ESTIMATE THE ABSOLUTE ERROR. 152C 153 RESG = 0.0E+00 154 FC = F(CENTR) 155 RESK = WGK(21)*FC 156 RESABS = ABS(RESK) 157 DO 10 J=1,10 158 JTW = J*2 159 ABSC = HLGTH*XGK(JTW) 160 FVAL1 = F(CENTR-ABSC) 161 FVAL2 = F(CENTR+ABSC) 162 FV1(JTW) = FVAL1 163 FV2(JTW) = FVAL2 164 FSUM = FVAL1+FVAL2 165 RESG = RESG+WG(J)*FSUM 166 RESK = RESK+WGK(JTW)*FSUM 167 RESABS = RESABS+WGK(JTW)*(ABS(FVAL1)+ABS(FVAL2)) 168 10 CONTINUE 169 DO 15 J = 1,10 170 JTWM1 = J*2-1 171 ABSC = HLGTH*XGK(JTWM1) 172 FVAL1 = F(CENTR-ABSC) 173 FVAL2 = F(CENTR+ABSC) 174 FV1(JTWM1) = FVAL1 175 FV2(JTWM1) = FVAL2 176 FSUM = FVAL1+FVAL2 177 RESK = RESK+WGK(JTWM1)*FSUM 178 RESABS = RESABS+WGK(JTWM1)*(ABS(FVAL1)+ABS(FVAL2)) 179 15 CONTINUE 180 RESKH = RESK*0.5E+00 181 RESASC = WGK(21)*ABS(FC-RESKH) 182 DO 20 J=1,20 183 RESASC = RESASC+WGK(J)*(ABS(FV1(J)-RESKH)+ABS(FV2(J)-RESKH)) 184 20 CONTINUE 185 RESULT = RESK*HLGTH 186 RESABS = RESABS*DHLGTH 187 RESASC = RESASC*DHLGTH 188 ABSERR = ABS((RESK-RESG)*HLGTH) 189 IF(RESASC.NE.0.0E+00.AND.ABSERR.NE.0.E+00) 190 1 ABSERR = RESASC*MIN(0.1E+01, 191 2 (0.2E+03*ABSERR/RESASC)**1.5E+00) 192 IF(RESABS.GT.UFLOW/(0.5E+02*EPMACH)) ABSERR = MAX 193 1 ((EPMACH*0.5E+02)*RESABS,ABSERR) 194 RETURN 195 END 196