estimate 1D integral on finite interval using a 31 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 modulus,
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 dqk31(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 31-point !! gauss-kronrod rule (resk), obtained by optimal !! addition of abscissae to the 15-point gauss !! rule (resg). real(wp), intent(out) :: Abserr !! estimate of the modulus of the modulus, !! 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 15-point gauss formula real(wp) :: resk !! result of the 31-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(15), fv2(15) 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(8), parameter :: wg = [ & 3.07532419961172683546283935772044177217e-2_wp, & 7.03660474881081247092674164506673384667e-2_wp, & 1.07159220467171935011869546685869303416e-1_wp, & 1.39570677926154314447804794511028322521e-1_wp, & 1.66269205816993933553200860481208811131e-1_wp, & 1.86161000015562211026800561866422824506e-1_wp, & 1.98431485327111576456118326443839324819e-1_wp, & 2.02578241925561272880620199967519314839e-1_wp] !! weights of the 15-point gauss rule real(wp), dimension(16), parameter :: xgk = [ & 9.98002298693397060285172840152271209073e-1_wp, & 9.87992518020485428489565718586612581147e-1_wp, & 9.67739075679139134257347978784337225283e-1_wp, & 9.37273392400705904307758947710209471244e-1_wp, & 8.97264532344081900882509656454495882832e-1_wp, & 8.48206583410427216200648320774216851366e-1_wp, & 7.90418501442465932967649294817947346862e-1_wp, & 7.24417731360170047416186054613938009631e-1_wp, & 6.50996741297416970533735895313274692547e-1_wp, & 5.70972172608538847537226737253910641238e-1_wp, & 4.85081863640239680693655740232350612866e-1_wp, & 3.94151347077563369897207370981045468363e-1_wp, & 2.99180007153168812166780024266388962662e-1_wp, & 2.01194093997434522300628303394596207813e-1_wp, & 1.01142066918717499027074231447392338787e-1_wp, & 0.00000000000000000000000000000000000000e0_wp] !! abscissae of the 31-point kronrod rule: !! !! * xgk(2), xgk(4), ... abscissae of the 15-point !! gauss rule !! * xgk(1), xgk(3), ... abscissae which are optimally !! added to the 15-point gauss rule real(wp), dimension(16), parameter :: wgk = [ & 5.37747987292334898779205143012764981831e-3_wp, & 1.50079473293161225383747630758072680946e-2_wp, & 2.54608473267153201868740010196533593973e-2_wp, & 3.53463607913758462220379484783600481226e-2_wp, & 4.45897513247648766082272993732796902233e-2_wp, & 5.34815246909280872653431472394302967716e-2_wp, & 6.20095678006706402851392309608029321904e-2_wp, & 6.98541213187282587095200770991474757860e-2_wp, & 7.68496807577203788944327774826590067221e-2_wp, & 8.30805028231330210382892472861037896016e-2_wp, & 8.85644430562117706472754436937743032123e-2_wp, & 9.31265981708253212254868727473457185619e-2_wp, & 9.66427269836236785051799076275893351367e-2_wp, & 9.91735987217919593323931734846031310596e-2_wp, & 1.00769845523875595044946662617569721916e-1_wp, & 1.01330007014791549017374792767492546771e-1_wp] !! weights of the 31-point kronrod rule centr = 0.5_wp*(a + b) hlgth = 0.5_wp*(b - a) dhlgth = abs(hlgth) ! compute the 31-point kronrod approximation to ! the integral, and estimate the absolute error. fc = f(centr) resg = wg(8)*fc resk = wgk(16)*fc Resabs = abs(resk) do j = 1, 7 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, 8 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(16)*abs(fc - reskh) do j = 1, 15 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 dqk31