subroutine dqk15w(f,w,p1,p2,p3,p4,kp,a,b,result,abserr, * resabs,resasc) c***begin prologue dqk15w c***date written 810101 (yymmdd) c***revision date 830518 (mmddyy) c***category no. h2a2a2 c***keywords 15-point gauss-kronrod rules c***author piessens,robert,appl. math. & progr. div. - k.u.leuven c de doncker,elise,appl. math. & progr. div. - k.u.leuven c***purpose to compute i = integral of f*w over (a,b), with error c estimate c j = integral of abs(f*w) over (a,b) c***description c c integration rules c standard fortran subroutine c double precision version c c parameters c on entry c f - double precision c function subprogram defining the integrand c function f(x). the actual name for f needs to be c declared e x t e r n a l in the driver program. c c w - double precision c function subprogram defining the integrand c weight function w(x). the actual name for w c needs to be declared e x t e r n a l in the c calling program. c c p1, p2, p3, p4 - double precision c parameters in the weight function c c kp - integer c key for indicating the type of weight function c c a - double precision c lower limit of integration c c b - double precision c upper limit of integration c c on return c result - double precision c approximation to the integral i c result is computed by applying the 15-point c kronrod rule (resk) obtained by optimal addition c of abscissae to the 7-point gauss rule (resg). c c abserr - double precision c estimate of the modulus of the absolute error, c which should equal or exceed abs(i-result) c c resabs - double precision c approximation to the integral of abs(f) c c resasc - double precision c approximation to the integral of abs(f-i/(b-a)) c c c***references (none) c***routines called d1mach c***end prologue dqk15w c double precision a,absc,absc1,absc2,abserr,b,centr,dabs,dhlgth, * dmax1,dmin1,d1mach,epmach,f,fc,fsum,fval1,fval2,fv1,fv2,hlgth, * p1,p2,p3,p4,resabs,resasc,resg,resk,reskh,result,uflow,w,wg,wgk, * xgk integer j,jtw,jtwm1,kp external f,w c dimension fv1(7),fv2(7),xgk(8),wgk(8),wg(4) c c the abscissae and weights are given for the interval (-1,1). c because of symmetry only the positive abscissae and their c corresponding weights are given. c c xgk - abscissae of the 15-point gauss-kronrod rule c xgk(2), xgk(4), ... abscissae of the 7-point c gauss rule c xgk(1), xgk(3), ... abscissae which are optimally c added to the 7-point gauss rule c c wgk - weights of the 15-point gauss-kronrod rule c c wg - weights of the 7-point gauss rule c data xgk(1),xgk(2),xgk(3),xgk(4),xgk(5),xgk(6),xgk(7),xgk(8)/ * 0.9914553711208126d+00, 0.9491079123427585d+00, * 0.8648644233597691d+00, 0.7415311855993944d+00, * 0.5860872354676911d+00, 0.4058451513773972d+00, * 0.2077849550078985d+00, 0.0000000000000000d+00/ c data wgk(1),wgk(2),wgk(3),wgk(4),wgk(5),wgk(6),wgk(7),wgk(8)/ * 0.2293532201052922d-01, 0.6309209262997855d-01, * 0.1047900103222502d+00, 0.1406532597155259d+00, * 0.1690047266392679d+00, 0.1903505780647854d+00, * 0.2044329400752989d+00, 0.2094821410847278d+00/ c data wg(1),wg(2),wg(3),wg(4)/ * 0.1294849661688697d+00, 0.2797053914892767d+00, * 0.3818300505051189d+00, 0.4179591836734694d+00/ c c c list of major variables c ----------------------- c c centr - mid point of the interval c hlgth - half-length of the interval c absc* - abscissa c fval* - function value c resg - result of the 7-point gauss formula c resk - result of the 15-point kronrod formula c reskh - approximation to the mean value of f*w over (a,b), c i.e. to i/(b-a) c c machine dependent constants c --------------------------- c c epmach is the largest relative spacing. c uflow is the smallest positive magnitude. c c***first executable statement dqk15w epmach = d1mach(4) uflow = d1mach(1) c centr = 0.5d+00*(a+b) hlgth = 0.5d+00*(b-a) dhlgth = dabs(hlgth) c c compute the 15-point kronrod approximation to the c integral, and estimate the error. c fc = f(centr)*w(centr,p1,p2,p3,p4,kp) resg = wg(4)*fc resk = wgk(8)*fc resabs = dabs(resk) do 10 j=1,3 jtw = j*2 absc = hlgth*xgk(jtw) absc1 = centr-absc absc2 = centr+absc fval1 = f(absc1)*w(absc1,p1,p2,p3,p4,kp) fval2 = f(absc2)*w(absc2,p1,p2,p3,p4,kp) fv1(jtw) = fval1 fv2(jtw) = fval2 fsum = fval1+fval2 resg = resg+wg(j)*fsum resk = resk+wgk(jtw)*fsum resabs = resabs+wgk(jtw)*(dabs(fval1)+dabs(fval2)) 10 continue do 15 j=1,4 jtwm1 = j*2-1 absc = hlgth*xgk(jtwm1) absc1 = centr-absc absc2 = centr+absc fval1 = f(absc1)*w(absc1,p1,p2,p3,p4,kp) fval2 = f(absc2)*w(absc2,p1,p2,p3,p4,kp) fv1(jtwm1) = fval1 fv2(jtwm1) = fval2 fsum = fval1+fval2 resk = resk+wgk(jtwm1)*fsum resabs = resabs+wgk(jtwm1)*(dabs(fval1)+dabs(fval2)) 15 continue reskh = resk*0.5d+00 resasc = wgk(8)*dabs(fc-reskh) do 20 j=1,7 resasc = resasc+wgk(j)*(dabs(fv1(j)-reskh)+dabs(fv2(j)-reskh)) 20 continue result = resk*hlgth resabs = resabs*dhlgth resasc = resasc*dhlgth abserr = dabs((resk-resg)*hlgth) if(resasc.ne.0.0d+00.and.abserr.ne.0.0d+00) * abserr = resasc*dmin1(0.1d+01,(0.2d+03*abserr/resasc)**1.5d+00) if(resabs.gt.uflow/(0.5d+02*epmach)) abserr = dmax1((epmach* * 0.5d+02)*resabs,abserr) return end