qk21.f

Go to the documentation of this file.
00001       subroutine qk21(f,a,b,result,abserr,resabs,resasc,ierr)
00002 c***begin prologue  qk21
00003 c***date written   800101   (yymmdd)
00004 c***revision date  830518   (yymmdd)
00005 c***category no.  h2a1a2
00006 c***keywords  21-point gauss-kronrod rules
00007 c***author  piessens,robert,appl. math. & progr. div. - k.u.leuven
00008 c           de doncker,elise,appl. math. & progr. div. - k.u.leuven
00009 c***purpose  to compute i = integral of f over (a,b), with error
00010 c                           estimate
00011 c                       j = integral of abs(f) over (a,b)
00012 c***description
00013 c
00014 c           integration rules
00015 c           standard fortran subroutine
00016 c           real version
00017 c
00018 c           parameters
00019 c            on entry
00020 c              f      - subroutine f(x,ierr,result) defining the integrand
00021 c                       function f(x). the actual name for f needs to be
00022 c                       declared e x t e r n a l in the driver program.
00023 c
00024 c              a      - real
00025 c                       lower limit of integration
00026 c
00027 c              b      - real
00028 c                       upper limit of integration
00029 c
00030 c            on return
00031 c              result - real
00032 c                       approximation to the integral i
00033 c                       result is computed by applying the 21-point
00034 c                       kronrod rule (resk) obtained by optimal addition
00035 c                       of abscissae to the 10-point gauss rule (resg).
00036 c
00037 c              abserr - real
00038 c                       estimate of the modulus of the absolute error,
00039 c                       which should not exceed abs(i-result)
00040 c
00041 c              resabs - real
00042 c                       approximation to the integral j
00043 c
00044 c              resasc - real
00045 c                       approximation to the integral of abs(f-i/(b-a))
00046 c                       over (a,b)
00047 c
00048 c***references  (none)
00049 c***routines called  r1mach
00050 c***end prologue  qk21
00051 c
00052       real a,absc,abserr,b,centr,dhlgth,epmach,fc,fsum,fval1,fval2,
00053      *  fv1,fv2,hlgth,resabs,resg,resk,reskh,result,r1mach,uflow,wg,wgk,
00054      *  xgk
00055       integer j,jtw,jtwm1
00056       external f
00057 c
00058       dimension fv1(10),fv2(10),wg(5),wgk(11),xgk(11)
00059 c
00060 c           the abscissae and weights are given for the interval (-1,1).
00061 c           because of symmetry only the positive abscissae and their
00062 c           corresponding weights are given.
00063 c
00064 c           xgk    - abscissae of the 21-point kronrod rule
00065 c                    xgk(2), xgk(4), ...  abscissae of the 10-point
00066 c                    gauss rule
00067 c                    xgk(1), xgk(3), ...  abscissae which are optimally
00068 c                    added to the 10-point gauss rule
00069 c
00070 c           wgk    - weights of the 21-point kronrod rule
00071 c
00072 c           wg     - weights of the 10-point gauss rule
00073 c
00074       data xgk(1),xgk(2),xgk(3),xgk(4),xgk(5),xgk(6),xgk(7),
00075      *  xgk(8),xgk(9),xgk(10),xgk(11)/
00076      *         0.9956571630258081e+00,     0.9739065285171717e+00,
00077      *     0.9301574913557082e+00,     0.8650633666889845e+00,
00078      *     0.7808177265864169e+00,     0.6794095682990244e+00,
00079      *     0.5627571346686047e+00,     0.4333953941292472e+00,
00080      *     0.2943928627014602e+00,     0.1488743389816312e+00,
00081      *     0.0000000000000000e+00/
00082 c
00083       data wgk(1),wgk(2),wgk(3),wgk(4),wgk(5),wgk(6),wgk(7),
00084      *  wgk(8),wgk(9),wgk(10),wgk(11)/
00085      *     0.1169463886737187e-01,     0.3255816230796473e-01,
00086      *     0.5475589657435200e-01,     0.7503967481091995e-01,
00087      *     0.9312545458369761e-01,     0.1093871588022976e+00,
00088      *     0.1234919762620659e+00,     0.1347092173114733e+00,
00089      *     0.1427759385770601e+00,     0.1477391049013385e+00,
00090      *     0.1494455540029169e+00/
00091 c
00092       data wg(1),wg(2),wg(3),wg(4),wg(5)/
00093      *     0.6667134430868814e-01,     0.1494513491505806e+00,
00094      *     0.2190863625159820e+00,     0.2692667193099964e+00,
00095      *     0.2955242247147529e+00/
00096 c
00097 c
00098 c           list of major variables
00099 c           -----------------------
00100 c
00101 c           centr  - mid point of the interval
00102 c           hlgth  - half-length of the interval
00103 c           absc   - abscissa
00104 c           fval*  - function value
00105 c           resg   - result of the 10-point gauss formula
00106 c           resk   - result of the 21-point kronrod formula
00107 c           reskh  - approximation to the mean value of f over (a,b),
00108 c                    i.e. to i/(b-a)
00109 c
00110 c
00111 c           machine dependent constants
00112 c           ---------------------------
00113 c
00114 c           epmach is the largest relative spacing.
00115 c           uflow is the smallest positive magnitude.
00116 c
00117 c***first executable statement  qk21
00118       epmach = r1mach(4)
00119       uflow = r1mach(1)
00120 c
00121       centr = 0.5e+00*(a+b)
00122       hlgth = 0.5e+00*(b-a)
00123       dhlgth = abs(hlgth)
00124 c
00125 c           compute the 21-point kronrod approximation to
00126 c           the integral, and estimate the absolute error.
00127 c
00128       resg = 0.0e+00
00129       call f(centr, ierr, fc)
00130       if (ierr .lt. 0) return
00131       resk = wgk(11)*fc
00132       resabs = abs(resk)
00133       do 10 j=1,5
00134         jtw = 2*j
00135         absc = hlgth*xgk(jtw)
00136         call f(centr-absc,ierr,fval1)
00137         if (ierr .lt. 0) return
00138         call f(centr+absc,ierr,fval2)
00139         if (ierr .lt. 0) return
00140         fv1(jtw) = fval1
00141         fv2(jtw) = fval2
00142         fsum = fval1+fval2
00143         resg = resg+wg(j)*fsum
00144         resk = resk+wgk(jtw)*fsum
00145         resabs = resabs+wgk(jtw)*(abs(fval1)+abs(fval2))
00146    10 continue
00147       do 15 j = 1,5
00148         jtwm1 = 2*j-1
00149         absc = hlgth*xgk(jtwm1)
00150         call f(centr-absc,ierr,fval1)
00151         if (ierr .lt. 0) return
00152         call f(centr+absc,ierr,fval2)
00153         if (ierr .lt. 0) return
00154         fv1(jtwm1) = fval1
00155         fv2(jtwm1) = fval2
00156         fsum = fval1+fval2
00157         resk = resk+wgk(jtwm1)*fsum
00158         resabs = resabs+wgk(jtwm1)*(abs(fval1)+abs(fval2))
00159    15 continue
00160       reskh = resk*0.5e+00
00161       resasc = wgk(11)*abs(fc-reskh)
00162       do 20 j=1,10
00163         resasc = resasc+wgk(j)*(abs(fv1(j)-reskh)+abs(fv2(j)-reskh))
00164    20 continue
00165       result = resk*hlgth
00166       resabs = resabs*dhlgth
00167       resasc = resasc*dhlgth
00168       abserr = abs((resk-resg)*hlgth)
00169       if(resasc.ne.0.0e+00.and.abserr.ne.0.0e+00)
00170      *  abserr = resasc*amin1(0.1e+01,
00171      *  (0.2e+03*abserr/resasc)**1.5e+00)
00172       if(resabs.gt.uflow/(0.5e+02*epmach)) abserr = amax1
00173      *  ((epmach*0.5e+02)*resabs,abserr)
00174       return
00175       end
 All Classes Files Functions Variables Typedefs Enumerations Enumerator Friends Defines