estimate 1D integral on finite interval using a 51 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 subroutine 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 not exceed |
||
| real(kind=wp), | intent(out) | :: | Resabs |
approximation to the integral j |
||
| real(kind=wp), | intent(out) | :: | Resasc |
approximation to the integral of |
subroutine dqk51(f, a, b, Result, Abserr, Resabs, Resasc) implicit none procedure(func) :: f !! function subroutine 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 51-point !! kronrod rule (resk) obtained by optimal addition !! of abscissae to the 25-point gauss rule (resg). real(wp), intent(out) :: Abserr !! estimate of the modulus of the absolute error, !! which should not 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))` !! over `(a,b)` 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 25-point gauss formula real(wp) :: resk !! result of the 51-point kronrod formula real(wp) :: reskh !! approximation to the mean value of `f` over `(a,b)`, i.e. to `i/(b-a)` real(wp) :: dhlgth, fc, fsum, fv1(25), fv2(25) integer :: j, jtw, jtwm1 ! 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(13), parameter :: wg = [ & 1.13937985010262879479029641132347736033e-2_wp, & 2.63549866150321372619018152952991449360e-2_wp, & 4.09391567013063126556234877116459536608e-2_wp, & 5.49046959758351919259368915404733241601e-2_wp, & 6.80383338123569172071871856567079685547e-2_wp, & 8.01407003350010180132349596691113022902e-2_wp, & 9.10282619829636498114972207028916533810e-2_wp, & 1.00535949067050644202206890392685826988e-1_wp, & 1.08519624474263653116093957050116619340e-1_wp, & 1.14858259145711648339325545869555808641e-1_wp, & 1.19455763535784772228178126512901047390e-1_wp, & 1.22242442990310041688959518945851505835e-1_wp, & 1.23176053726715451203902873079050142438e-1_wp] !! weights of the 25-point gauss rule real(wp), dimension(26), parameter :: xgk = [ & 9.99262104992609834193457486540340593705e-1_wp, & 9.95556969790498097908784946893901617258e-1_wp, & 9.88035794534077247637331014577406227072e-1_wp, & 9.76663921459517511498315386479594067745e-1_wp, & 9.61614986425842512418130033660167241692e-1_wp, & 9.42974571228974339414011169658470531905e-1_wp, & 9.20747115281701561746346084546330631575e-1_wp, & 8.94991997878275368851042006782804954175e-1_wp, & 8.65847065293275595448996969588340088203e-1_wp, & 8.33442628760834001421021108693569569461e-1_wp, & 7.97873797998500059410410904994306569409e-1_wp, & 7.59259263037357630577282865204360976388e-1_wp, & 7.17766406813084388186654079773297780598e-1_wp, & 6.73566368473468364485120633247622175883e-1_wp, & 6.26810099010317412788122681624517881020e-1_wp, & 5.77662930241222967723689841612654067396e-1_wp, & 5.26325284334719182599623778158010178037e-1_wp, & 4.73002731445714960522182115009192041332e-1_wp, & 4.17885382193037748851814394594572487093e-1_wp, & 3.61172305809387837735821730127640667422e-1_wp, & 3.03089538931107830167478909980339329200e-1_wp, & 2.43866883720988432045190362797451586406e-1_wp, & 1.83718939421048892015969888759528415785e-1_wp, & 1.22864692610710396387359818808036805532e-1_wp, & 6.15444830056850788865463923667966312817e-2_wp, & 0.00000000000000000000000000000000000000e0_wp] !! abscissae of the 51-point kronrod rule !! !! * xgk(2), xgk(4), ... abscissae of the 25-point !! gauss rule !! * xgk(1), xgk(3), ... abscissae which are optimally !! added to the 25-point gauss rule real(wp), dimension(26), parameter :: wgk = [ & 1.98738389233031592650785188284340988943e-3_wp, & 5.56193213535671375804023690106552207018e-3_wp, & 9.47397338617415160720771052365532387165e-3_wp, & 1.32362291955716748136564058469762380776e-2_wp, & 1.68478177091282982315166675363363158404e-2_wp, & 2.04353711458828354565682922359389736788e-2_wp, & 2.40099456069532162200924891648810813929e-2_wp, & 2.74753175878517378029484555178110786148e-2_wp, & 3.07923001673874888911090202152285856009e-2_wp, & 3.40021302743293378367487952295512032257e-2_wp, & 3.71162714834155435603306253676198759960e-2_wp, & 4.00838255040323820748392844670756464014e-2_wp, & 4.28728450201700494768957924394951611020e-2_wp, & 4.55029130499217889098705847526603930437e-2_wp, & 4.79825371388367139063922557569147549836e-2_wp, & 5.02776790807156719633252594334400844406e-2_wp, & 5.23628858064074758643667121378727148874e-2_wp, & 5.42511298885454901445433704598756068261e-2_wp, & 5.59508112204123173082406863827473468203e-2_wp, & 5.74371163615678328535826939395064719948e-2_wp, & 5.86896800223942079619741758567877641398e-2_wp, & 5.97203403241740599790992919325618538354e-2_wp, & 6.05394553760458629453602675175654271623e-2_wp, & 6.11285097170530483058590304162927119227e-2_wp, & 6.14711898714253166615441319652641775865e-2_wp, & 6.15808180678329350787598242400645531904e-2_wp] !! weights of the 51-point kronrod rule. centr = 0.5_wp*(a + b) hlgth = 0.5_wp*(b - a) dhlgth = abs(hlgth) ! compute the 51-point kronrod approximation to ! the integral, and estimate the absolute error. fc = f(centr) resg = wg(13)*fc resk = wgk(26)*fc Resabs = abs(resk) do j = 1, 12 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, 13 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(26)*abs(fc - reskh) do j = 1, 25 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 dqk51