estimate 1D integral on finite interval using a 21 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 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 dqk21(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 21-point !! kronrod rule (resk) obtained by optimal addition !! of abscissae to the 10-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) :: dhlgth, fc, fsum, fv1(10), fv2(10) 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 10-point gauss formula real(wp) :: resk !! result of the 21-point kronrod formula 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(5), parameter :: wg = [ & 6.66713443086881375935688098933317928579e-2_wp, & 1.49451349150580593145776339657697332403e-1_wp, & 2.19086362515982043995534934228163192459e-1_wp, & 2.69266719309996355091226921569469352860e-1_wp, & 2.95524224714752870173892994651338329421e-1_wp] !! weights of the 10-point gauss rule real(wp), dimension(11), parameter :: xgk = [ & 9.95657163025808080735527280689002847921e-1_wp, & 9.73906528517171720077964012084452053428e-1_wp, & 9.30157491355708226001207180059508346225e-1_wp, & 8.65063366688984510732096688423493048528e-1_wp, & 7.80817726586416897063717578345042377163e-1_wp, & 6.79409568299024406234327365114873575769e-1_wp, & 5.62757134668604683339000099272694140843e-1_wp, & 4.33395394129247190799265943165784162200e-1_wp, & 2.94392862701460198131126603103865566163e-1_wp, & 1.48874338981631210884826001129719984618e-1_wp, & 0.00000000000000000000000000000000000000e0_wp] !! abscissae of the 21-point kronrod rule: !! !! * xgk(2), xgk(4), ... abscissae of the 10-point !! gauss rule !! * xgk(1), xgk(3), ... abscissae which are optimally !! added to the 10-point gauss rule real(wp), dimension(11), parameter :: wgk = [ & 1.16946388673718742780643960621920483962e-2_wp, & 3.25581623079647274788189724593897606174e-2_wp, & 5.47558965743519960313813002445801763737e-2_wp, & 7.50396748109199527670431409161900093952e-2_wp, & 9.31254545836976055350654650833663443900e-2_wp, & 1.09387158802297641899210590325804960272e-1_wp, & 1.23491976262065851077958109831074159512e-1_wp, & 1.34709217311473325928054001771706832761e-1_wp, & 1.42775938577060080797094273138717060886e-1_wp, & 1.47739104901338491374841515972068045524e-1_wp, & 1.49445554002916905664936468389821203745e-1_wp] !! weights of the 21-point kronrod rule centr = 0.5_wp*(a + b) hlgth = 0.5_wp*(b - a) dhlgth = abs(hlgth) ! compute the 21-point kronrod approximation to ! the integral, and estimate the absolute error. resg = 0.0_wp fc = f(centr) resk = wgk(11)*fc Resabs = abs(resk) do j = 1, 5 jtw = 2*j 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, 5 jtwm1 = 2*j - 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(11)*abs(fc - reskh) do j = 1, 10 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 dqk21