estimate 1D integral on finite interval using a 61 point gauss-kronrod rule and give error estimate, non-automatic
to compute i = integral of f over (a,b) with error
estimate j = integral of abs(f) over (a,b).
| Type | Intent | Optional | Attributes | Name | ||
|---|---|---|---|---|---|---|
| procedure(func) | :: | f |
function subprogram defining the integrand function |
|||
| real(kind=wp), | intent(in) | :: | a |
lower limit of integration |
||
| real(kind=wp), | intent(in) | :: | b |
upper limit of integration |
||
| real(kind=wp), | intent(out) | :: | Result |
approximation to the integral i
|
||
| real(kind=wp), | intent(out) | :: | Abserr |
estimate of the modulus of the absolute error,
which should equal or exceed |
||
| real(kind=wp), | intent(out) | :: | Resabs |
approximation to the integral j |
||
| real(kind=wp), | intent(out) | :: | Resasc |
approximation to the integral of |
subroutine dqk61(f, a, b, Result, Abserr, Resabs, Resasc) implicit none procedure(func) :: f !! function subprogram defining the integrand function `f(x)`. real(wp), intent(in) :: a !! lower limit of integration real(wp), intent(in) :: b !! upper limit of integration real(wp), intent(out) :: Result !! approximation to the integral i !! `result` is computed by applying the 61-point !! kronrod rule (resk) obtained by optimal addition of !! abscissae to the 30-point gauss rule (resg). real(wp), intent(out) :: Abserr !! estimate of the modulus of the absolute error, !! which should equal or exceed `abs(i-result)` real(wp), intent(out) :: Resabs !! approximation to the integral j real(wp), intent(out) :: Resasc !! approximation to the integral of `abs(f-i/(b-a))` real(wp) :: dhlgth, fc, fsum, fv1(30), fv2(30) integer :: j, jtw, jtwm1 real(wp) :: centr !! mid point of the interval real(wp) :: hlgth !! half-length of the interval real(wp) :: absc !! abscissa real(wp) :: fval1 !! function value real(wp) :: fval2 !! function value real(wp) :: resg !! result of the 30-point gauss rule real(wp) :: resk !! result of the 61-point kronrod rule real(wp) :: reskh !! approximation to the mean value of `f` over `(a,b)`, i.e. to `i/(b-a)` ! the abscissae and weights are given for the ! interval (-1,1). because of symmetry only the positive ! abscissae and their corresponding weights are given. real(wp), dimension(15), parameter :: wg = [ & 7.96819249616660561546588347467362245048e-3_wp, & 1.84664683110909591423021319120472690962e-2_wp, & 2.87847078833233693497191796112920436396e-2_wp, & 3.87991925696270495968019364463476920332e-2_wp, & 4.84026728305940529029381404228075178153e-2_wp, & 5.74931562176190664817216894020561287971e-2_wp, & 6.59742298821804951281285151159623612374e-2_wp, & 7.37559747377052062682438500221907341538e-2_wp, & 8.07558952294202153546949384605297308759e-2_wp, & 8.68997872010829798023875307151257025768e-2_wp, & 9.21225222377861287176327070876187671969e-2_wp, & 9.63687371746442596394686263518098650964e-2_wp, & 9.95934205867952670627802821035694765299e-2_wp, & 1.01762389748405504596428952168554044633e-1_wp, & 1.02852652893558840341285636705415043868e-1_wp] !! weights of the 30-point gauss rule real(wp), dimension(31), parameter :: xgk = [ & 9.99484410050490637571325895705810819469e-1_wp, & 9.96893484074649540271630050918695283341e-1_wp, & 9.91630996870404594858628366109485724851e-1_wp, & 9.83668123279747209970032581605662801940e-1_wp, & 9.73116322501126268374693868423706884888e-1_wp, & 9.60021864968307512216871025581797662930e-1_wp, & 9.44374444748559979415831324037439121586e-1_wp, & 9.26200047429274325879324277080474004086e-1_wp, & 9.05573307699907798546522558925958319569e-1_wp, & 8.82560535792052681543116462530225590057e-1_wp, & 8.57205233546061098958658510658943856821e-1_wp, & 8.29565762382768397442898119732501916439e-1_wp, & 7.99727835821839083013668942322683240736e-1_wp, & 7.67777432104826194917977340974503131695e-1_wp, & 7.33790062453226804726171131369527645669e-1_wp, & 6.97850494793315796932292388026640068382e-1_wp, & 6.60061064126626961370053668149270753038e-1_wp, & 6.20526182989242861140477556431189299207e-1_wp, & 5.79345235826361691756024932172540495907e-1_wp, & 5.36624148142019899264169793311072794164e-1_wp, & 4.92480467861778574993693061207708795644e-1_wp, & 4.47033769538089176780609900322854000162e-1_wp, & 4.00401254830394392535476211542660633611e-1_wp, & 3.52704725530878113471037207089373860654e-1_wp, & 3.04073202273625077372677107199256553531e-1_wp, & 2.54636926167889846439805129817805107883e-1_wp, & 2.04525116682309891438957671002024709524e-1_wp, & 1.53869913608583546963794672743255920419e-1_wp, & 1.02806937966737030147096751318000592472e-1_wp, & 5.14718425553176958330252131667225737491e-2_wp, & 0.00000000000000000000000000000000000000e0_wp] !! abscissae of the 61-point kronrod rule: !! !! * `xgk(2), xgk(4)` ... abscissae of the 30-point !! gauss rule !! * `xgk(1), xgk(3)` ... optimally added abscissae !! to the 30-point gauss rule real(wp), dimension(31), parameter :: wgk = [ & 1.38901369867700762455159122675969968105e-3, & 3.89046112709988405126720184451550327852e-3, & 6.63070391593129217331982636975016813363e-3, & 9.27327965951776342844114689202436042127e-3, & 1.18230152534963417422328988532505928963e-2, & 1.43697295070458048124514324435800101958e-2, & 1.69208891890532726275722894203220923686e-2, & 1.94141411939423811734089510501284558514e-2, & 2.18280358216091922971674857383389934015e-2, & 2.41911620780806013656863707252320267604e-2, & 2.65099548823331016106017093350754143665e-2, & 2.87540487650412928439787853543342111447e-2, & 3.09072575623877624728842529430922726353e-2, & 3.29814470574837260318141910168539275106e-2, & 3.49793380280600241374996707314678750972e-2, & 3.68823646518212292239110656171359677370e-2, & 3.86789456247275929503486515322810502509e-2, & 4.03745389515359591119952797524681142161e-2, & 4.19698102151642461471475412859697577901e-2, & 4.34525397013560693168317281170732580746e-2, & 4.48148001331626631923555516167232437574e-2, & 4.60592382710069881162717355593735805947e-2, & 4.71855465692991539452614781810994864829e-2, & 4.81858617570871291407794922983045926058e-2, & 4.90554345550297788875281653672381736059e-2, & 4.97956834270742063578115693799423285392e-2, & 5.04059214027823468408930856535850289022e-2, & 5.08817958987496064922974730498046918534e-2, & 5.12215478492587721706562826049442082511e-2, & 5.14261285374590259338628792157812598296e-2, & 5.14947294294515675583404336470993075327e-2] !! weights of the 61-point kronrod rule centr = 0.5_wp*(b + a) hlgth = 0.5_wp*(b - a) dhlgth = abs(hlgth) ! compute the 61-point kronrod approximation to the ! integral, and estimate the absolute error. resg = 0.0_wp fc = f(centr) resk = wgk(31)*fc Resabs = abs(resk) do j = 1, 15 jtw = j*2 absc = hlgth*xgk(jtw) fval1 = f(centr - absc) fval2 = f(centr + absc) fv1(jtw) = fval1 fv2(jtw) = fval2 fsum = fval1 + fval2 resg = resg + wg(j)*fsum resk = resk + wgk(jtw)*fsum Resabs = Resabs + wgk(jtw)*(abs(fval1) + abs(fval2)) end do do j = 1, 15 jtwm1 = j*2 - 1 absc = hlgth*xgk(jtwm1) fval1 = f(centr - absc) fval2 = f(centr + absc) fv1(jtwm1) = fval1 fv2(jtwm1) = fval2 fsum = fval1 + fval2 resk = resk + wgk(jtwm1)*fsum Resabs = Resabs + wgk(jtwm1)*(abs(fval1) + abs(fval2)) end do reskh = resk*0.5_wp Resasc = wgk(31)*abs(fc - reskh) do j = 1, 30 Resasc = Resasc + wgk(j) & *(abs(fv1(j) - reskh) + abs(fv2(j) - reskh)) end do Result = resk*hlgth Resabs = Resabs*dhlgth Resasc = Resasc*dhlgth Abserr = abs((resk - resg)*hlgth) if (Resasc /= 0.0_wp .and. Abserr /= 0.0_wp) & Abserr = Resasc*min(1.0_wp, (200.0_wp*Abserr/Resasc)**1.5_wp) if (Resabs > uflow/(50.0_wp*epmach)) & Abserr = max((epmach*50.0_wp)*Resabs, Abserr) end subroutine dqk61