subroutine qc25s(f,a,b,bl,br,alfa,beta,ri,rj,rg,rh,result, * abserr,resasc,integr,nev) c***begin prologue qc25s c***date written 810101 (yymmdd) c***revision date 830518 (yymmdd) c***category no. h2a2a2 c***keywords 25-point clenshaw-curtis integration 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 (bl,br), with error c estimate, where the weight function w has a singular c behaviour of algebraico-logarithmic type at the points c a and/or b. (bl,br) is a part of (a,b). c***description c c integration rules for integrands having algebraico-logarithmic c end point singularities c standard fortran subroutine c real version c c parameters c f - real c function subprogram defining the integrand c f(x). the actual name for f needs to be declared c e x t e r n a l in the driver program. c c a - real c left end point of the original interval c c b - real c right end point of the original interval, b.gt.a c c bl - real c lower limit of integration, bl.ge.a c c br - real c upper limit of integration, br.le.b c c alfa - real c parameter in the weight function c c beta - real c parameter in the weight function c c ri,rj,rg,rh - real c modified chebyshev moments for the application c of the generalized clenshaw-curtis c method (computed in subroutine dqmomo) c c result - real c approximation to the integral c result is computed by using a generalized c clenshaw-curtis method if b1 = a or br = b. c in all other cases the 15-point kronrod c rule is applied, obtained by optimal addition of c abscissae to the 7-point gauss rule. c c abserr - real c estimate of the modulus of the absolute error, c which should equal or exceed abs(i-result) c c resasc - real c approximation to the integral of abs(f*w-i/(b-a)) c c integr - integer c which determines the weight function c = 1 w(x) = (x-a)**alfa*(b-x)**beta c = 2 w(x) = (x-a)**alfa*(b-x)**beta*log(x-a) c = 3 w(x) = (x-a)**alfa*(b-x)**beta*log(b-x) c = 4 w(x) = (x-a)**alfa*(b-x)**beta*log(x-a)* c log(b-x) c c nev - integer c number of integrand evaluations c c***references (none) c***routines called qcheb,qk15w c***end prologue qc25s c real a,abserr,alfa,b,beta,bl,br,centr,cheb12,cheb24, * dc,f,factor,fix,fval,hlgth,resabs,resasc, * result,res12,res24,rg,rh,ri,rj,u,qwgts,x integer i,integr,isym,nev c dimension cheb12(13),cheb24(25),fval(25),rg(25),rh(25),ri(25), * rj(25),x(11) c external f,qwgts c c the vector x contains the values cos(k*pi/24) c k = 1, ..., 11, to be used for the computation of the c chebyshev series expansion of f. c data x(1),x(2),x(3),x(4),x(5),x(6),x(7),x(8),x(9),x(10), * x(11)/ * 0.9914448613738104e+00, 0.9659258262890683e+00, * 0.9238795325112868e+00, 0.8660254037844386e+00, * 0.7933533402912352e+00, 0.7071067811865475e+00, * 0.6087614290087206e+00, 0.5000000000000000e+00, * 0.3826834323650898e+00, 0.2588190451025208e+00, * 0.1305261922200516e+00/ c c list of major variables c ----------------------- c c fval - value of the function f at the points c (br-bl)*0.5*cos(k*pi/24)+(br+bl)*0.5 c k = 0, ..., 24 c cheb12 - coefficients of the chebyshev series expansion c of degree 12, for the function f, in the c interval (bl,br) c cheb24 - coefficients of the chebyshev series expansion c of degree 24, for the function f, in the c interval (bl,br) c res12 - approximation to the integral obtained from cheb12 c res24 - approximation to the integral obtained from cheb24 c qwgts - external function subprogram defining c the four possible weight functions c hlgth - half-length of the interval (bl,br) c centr - mid point of the interval (bl,br) c c***first executable statement qc25s nev = 25 if(bl.eq.a.and.(alfa.ne.0.0e+00.or.integr.eq.2.or.integr.eq.4)) * go to 10 if(br.eq.b.and.(beta.ne.0.0e+00.or.integr.eq.3.or.integr.eq.4)) * go to 140 c c if a.gt.bl and b.lt.br, apply the 15-point gauss-kronrod c scheme. c c call qk15w(f,qwgts,a,b,alfa,beta,integr,bl,br, * result,abserr,resabs,resasc) nev = 15 go to 270 c c this part of the program is executed only if a = bl. c ---------------------------------------------------- c c compute the chebyshev series expansion of the c following function c f1 = (0.5*(b+b-br-a)-0.5*(br-a)*x)**beta c *f(0.5*(br-a)*x+0.5*(br+a)) c 10 hlgth = 0.5e+00*(br-bl) centr = 0.5e+00*(br+bl) fix = b-centr fval(1) = 0.5e+00*f(hlgth+centr)*(fix-hlgth)**beta fval(13) = f(centr)*(fix**beta) fval(25) = 0.5e+00*f(centr-hlgth)*(fix+hlgth)**beta do 20 i=2,12 u = hlgth*x(i-1) isym = 26-i fval(i) = f(u+centr)*(fix-u)**beta fval(isym) = f(centr-u)*(fix+u)**beta 20 continue factor = hlgth**(alfa+0.1e+01) result = 0.0e+00 abserr = 0.0e+00 res12 = 0.0e+00 res24 = 0.0e+00 if(integr.gt.2) go to 70 call qcheb(x,fval,cheb12,cheb24) c c integr = 1 (or 2) c do 30 i=1,13 res12 = res12+cheb12(i)*ri(i) res24 = res24+cheb24(i)*ri(i) 30 continue do 40 i=14,25 res24 = res24+cheb24(i)*ri(i) 40 continue if(integr.eq.1) go to 130 c c integr = 2 c dc = alog(br-bl) result = res24*dc abserr = abs((res24-res12)*dc) res12 = 0.0e+00 res24 = 0.0e+00 do 50 i=1,13 res12 = res12+cheb12(i)*rg(i) res24 = res12+cheb24(i)*rg(i) 50 continue do 60 i=14,25 res24 = res24+cheb24(i)*rg(i) 60 continue go to 130 c c compute the chebyshev series expansion of the c following function c f4 = f1*log(0.5*(b+b-br-a)-0.5*(br-a)*x) c 70 fval(1) = fval(1)*alog(fix-hlgth) fval(13) = fval(13)*alog(fix) fval(25) = fval(25)*alog(fix+hlgth) do 80 i=2,12 u = hlgth*x(i-1) isym = 26-i fval(i) = fval(i)*alog(fix-u) fval(isym) = fval(isym)*alog(fix+u) 80 continue call qcheb(x,fval,cheb12,cheb24) c c integr = 3 (or 4) c do 90 i=1,13 res12 = res12+cheb12(i)*ri(i) res24 = res24+cheb24(i)*ri(i) 90 continue do 100 i=14,25 res24 = res24+cheb24(i)*ri(i) 100 continue if(integr.eq.3) go to 130 c c integr = 4 c dc = alog(br-bl) result = res24*dc abserr = abs((res24-res12)*dc) res12 = 0.0e+00 res24 = 0.0e+00 do 110 i=1,13 res12 = res12+cheb12(i)*rg(i) res24 = res24+cheb24(i)*rg(i) 110 continue do 120 i=14,25 res24 = res24+cheb24(i)*rg(i) 120 continue 130 result = (result+res24)*factor abserr = (abserr+abs(res24-res12))*factor go to 270 c c this part of the program is executed only if b = br. c ---------------------------------------------------- c c compute the chebyshev series expansion of the c following function c f2 = (0.5*(b+bl-a-a)+0.5*(b-bl)*x)**alfa c *f(0.5*(b-bl)*x+0.5*(b+bl)) c 140 hlgth = 0.5e+00*(br-bl) centr = 0.5e+00*(br+bl) fix = centr-a fval(1) = 0.5e+00*f(hlgth+centr)*(fix+hlgth)**alfa fval(13) = f(centr)*(fix**alfa) fval(25) = 0.5e+00*f(centr-hlgth)*(fix-hlgth)**alfa do 150 i=2,12 u = hlgth*x(i-1) isym = 26-i fval(i) = f(u+centr)*(fix+u)**alfa fval(isym) = f(centr-u)*(fix-u)**alfa 150 continue factor = hlgth**(beta+0.1e+01) result = 0.0e+00 abserr = 0.0e+00 res12 = 0.0e+00 res24 = 0.0e+00 if(integr.eq.2.or.integr.eq.4) go to 200 c c integr = 1 (or 3) c call qcheb(x,fval,cheb12,cheb24) do 160 i=1,13 res12 = res12+cheb12(i)*rj(i) res24 = res24+cheb24(i)*rj(i) 160 continue do 170 i=14,25 res24 = res24+cheb24(i)*rj(i) 170 continue if(integr.eq.1) go to 260 c c integr = 3 c dc = alog(br-bl) result = res24*dc abserr = abs((res24-res12)*dc) res12 = 0.0e+00 res24 = 0.0e+00 do 180 i=1,13 res12 = res12+cheb12(i)*rh(i) res24 = res24+cheb24(i)*rh(i) 180 continue do 190 i=14,25 res24 = res24+cheb24(i)*rh(i) 190 continue go to 260 c c compute the chebyshev series expansion of the c following function c f3 = f2*log(0.5*(b-bl)*x+0.5*(b+bl-a-a)) c 200 fval(1) = fval(1)*alog(hlgth+fix) fval(13) = fval(13)*alog(fix) fval(25) = fval(25)*alog(fix-hlgth) do 210 i=2,12 u = hlgth*x(i-1) isym = 26-i fval(i) = fval(i)*alog(u+fix) fval(isym) = fval(isym)*alog(fix-u) 210 continue call qcheb(x,fval,cheb12,cheb24) c c integr = 2 (or 4) c do 220 i=1,13 res12 = res12+cheb12(i)*rj(i) res24 = res24+cheb24(i)*rj(i) 220 continue do 230 i=14,25 res24 = res24+cheb24(i)*rj(i) 230 continue if(integr.eq.2) go to 260 dc = alog(br-bl) result = res24*dc abserr = abs((res24-res12)*dc) res12 = 0.0e+00 res24 = 0.0e+00 c c integr = 4 c do 240 i=1,13 res12 = res12+cheb12(i)*rh(i) res24 = res24+cheb24(i)*rh(i) 240 continue do 250 i=14,25 res24 = res24+cheb24(i)*rh(i) 250 continue 260 result = (result+res24)*factor abserr = (abserr+abs(res24-res12))*factor 270 return end