estimate 1D integral on (semi)infinite interval using a 15 point gauss-kronrod quadrature rule, non-automatic
the original (infinite integration range is mapped onto the interval (0,1) and (a,b) is a part of (0,1). it is the purpose to compute:
(a,b),(a,b).| Type | Intent | Optional | Attributes | Name | ||
|---|---|---|---|---|---|---|
| procedure(func) | :: | f |
function subprogram defining the integrand function |
|||
| real(kind=wp), | intent(in) | :: | Boun |
finite bound of original integration range (set to zero if inf = +2) |
||
| integer, | intent(in) | :: | Inf |
the integral is computed as the sum of two
integrals, one over |
||
| real(kind=wp), | intent(in) | :: | a |
lower limit for integration over subrange of |
||
| real(kind=wp), | intent(in) | :: | b |
upper limit for integration over subrange of |
||
| 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 dqk15i(f, Boun, Inf, a, b, Result, Abserr, Resabs, Resasc) implicit none procedure(func) :: f !! function subprogram defining the integrand function `f(x)`. real(wp), intent(in) :: Boun !! finite bound of original integration !! range (set to zero if inf = +2) real(wp), intent(in) :: a !! lower limit for integration over subrange of `(0,1)` real(wp), intent(in) :: b !! upper limit for integration over subrange of `(0,1)` integer, intent(in) :: Inf !! * if inf = -1, the original interval is !! `(-infinity,bound)`, !! * if inf = +1, the original interval is !! `(bound,+infinity)`, !! * if inf = +2, the original interval is !! `(-infinity,+infinity)` and !! !! the integral is computed as the sum of two !! integrals, one over `(-infinity,0)` and one over !! `(0,+infinity)`. real(wp), intent(out) :: Result !! approximation to the integral i. !! `result` is computed by applying the 15-point !! kronrod rule(resk) obtained by optimal addition !! of abscissae to the 7-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((transformed integrand)-i/(b-a))` over `(a,b)` real(wp) :: absc, dinf, fc, fsum, fv1(7), fv2(7) integer :: j real(wp) :: centr !! mid point of the interval real(wp) :: hlgth !! half-length of the interval real(wp) :: absc1 !! abscissa real(wp) :: absc2 !! abscissa real(wp) :: tabsc1 !! transformed abscissa real(wp) :: tabsc2 !! transformed abscissa real(wp) :: fval1 !! function value real(wp) :: fval2 !! function value real(wp) :: resg !! result of the 7-point gauss formula real(wp) :: resk !! result of the 15-point kronrod formula real(wp) :: reskh !! approximation to the mean value of the transformed !! integrand over `(a,b)`, i.e. to `i/(b-a)` ! the abscissae and weights are supplied for the interval ! (-1,1). because of symmetry only the positive abscissae and ! their corresponding weights are given. real(wp), dimension(8), parameter :: wg = [ & 0.00000000000000000000000000000000000000e0_wp, & 1.29484966168869693270611432679082018329e-1_wp, & 0.00000000000000000000000000000000000000e0_wp, & 2.79705391489276667901467771423779582487e-1_wp, & 0.00000000000000000000000000000000000000e0_wp, & 3.81830050505118944950369775488975133878e-1_wp, & 0.00000000000000000000000000000000000000e0_wp, & 4.17959183673469387755102040816326530612e-1_wp] !! weights of the 7-point gauss rule, corresponding !! to the abscissae `xgk(2), xgk(4), ...`. !! `wg(1), wg(3), ...` are set to zero. real(wp), dimension(8), parameter :: xgk = [ & 9.91455371120812639206854697526328516642e-1_wp, & 9.49107912342758524526189684047851262401e-1_wp, & 8.64864423359769072789712788640926201211e-1_wp, & 7.41531185599394439863864773280788407074e-1_wp, & 5.86087235467691130294144838258729598437e-1_wp, & 4.05845151377397166906606412076961463347e-1_wp, & 2.07784955007898467600689403773244913480e-1_wp, & 0.00000000000000000000000000000000000000e0_wp] !! abscissae of the 15-point kronrod rule: !! !! * xgk(2), xgk(4), ... abscissae of the 7-point !! gauss rule !! * xgk(1), xgk(3), ... abscissae which are optimally !! added to the 7-point gauss rule real(wp), dimension(8), parameter :: wgk = [ & 2.29353220105292249637320080589695919936e-2_wp, & 6.30920926299785532907006631892042866651e-2_wp, & 1.04790010322250183839876322541518017444e-1_wp, & 1.40653259715525918745189590510237920400e-1_wp, & 1.69004726639267902826583426598550284106e-1_wp, & 1.90350578064785409913256402421013682826e-1_wp, & 2.04432940075298892414161999234649084717e-1_wp, & 2.09482141084727828012999174891714263698e-1_wp] !! weights of the 15-point kronrod rule dinf = min(1, Inf) centr = 0.5_wp*(a + b) hlgth = 0.5_wp*(b - a) tabsc1 = Boun + dinf*(1.0_wp - centr)/centr fval1 = f(tabsc1) if (Inf == 2) fval1 = fval1 + f(-tabsc1) fc = (fval1/centr)/centr ! compute the 15-point kronrod approximation to ! the integral, and estimate the error. resg = wg(8)*fc resk = wgk(8)*fc Resabs = abs(resk) do j = 1, 7 absc = hlgth*xgk(j) absc1 = centr - absc absc2 = centr + absc tabsc1 = Boun + dinf*(1.0_wp - absc1)/absc1 tabsc2 = Boun + dinf*(1.0_wp - absc2)/absc2 fval1 = f(tabsc1) fval2 = f(tabsc2) if (Inf == 2) then fval1 = fval1 + f(-tabsc1) fval2 = fval2 + f(-tabsc2) end if fval1 = (fval1/absc1)/absc1 fval2 = (fval2/absc2)/absc2 fv1(j) = fval1 fv2(j) = fval2 fsum = fval1 + fval2 resg = resg + wg(j)*fsum resk = resk + wgk(j)*fsum Resabs = Resabs + wgk(j)*(abs(fval1) + abs(fval2)) end do reskh = resk*0.5_wp Resasc = wgk(8)*abs(fc - reskh) do j = 1, 7 Resasc = Resasc + wgk(j) & *(abs(fv1(j) - reskh) + abs(fv2(j) - reskh)) end do Result = resk*hlgth Resasc = Resasc*hlgth Resabs = Resabs*hlgth 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 dqk15i