dqc25s Subroutine

public subroutine dqc25s(f, a, b, Bl, Br, Alfa, Beta, Ri, Rj, Rg, Rh, Result, Abserr, Resasc, Integr, Nev)

25-point clenshaw-curtis integration

to compute i = integral of f*w over (bl,br), with error estimate, where the weight function w has a singular behaviour of algebraico-logarithmic type at the points a and/or b. (bl,br) is a part of (a,b).

History

  • QUADPACK: date written 810101, revision date 830518 (yymmdd)

Arguments

Type IntentOptional Attributes Name
procedure(func) :: f

function subprogram defining the integrand f(x).
real(kind=wp), intent(in) :: a

left end point of the original interval
real(kind=wp), intent(in) :: b

right end point of the original interval, b>a
real(kind=wp), intent(in) :: Bl

lower limit of integration, bl>=a
real(kind=wp), intent(in) :: Br

upper limit of integration, br<=b
real(kind=wp), intent(in) :: Alfa

parameter in the weight function
real(kind=wp), intent(in) :: Beta

parameter in the weight function
real(kind=wp), intent(in) :: Ri(25)

modified chebyshev moments for the application of the generalized clenshaw-curtis method (computed in subroutine dqmomo)
real(kind=wp), intent(in) :: Rj(25)

modified chebyshev moments for the application of the generalized clenshaw-curtis method (computed in subroutine dqmomo)
real(kind=wp), intent(in) :: Rg(25)

modified chebyshev moments for the application of the generalized clenshaw-curtis method (computed in subroutine dqmomo)
real(kind=wp), intent(in) :: Rh(25)

modified chebyshev moments for the application of the generalized clenshaw-curtis method (computed in subroutine dqmomo)
real(kind=wp), intent(out) :: Result

approximation to the integral result is computed by using a generalized clenshaw-curtis method if b1 = a or br = b. in all other cases the 15-point kronrod rule is applied, obtained by optimal addition of abscissae to the 7-point gauss rule.
real(kind=wp), intent(out) :: Abserr

estimate of the modulus of the absolute error, which should equal or exceed abs(i-result)
real(kind=wp), intent(out) :: Resasc

approximation to the integral of abs(f*w-i/(b-a))
integer, intent(in) :: Integr

which determines the weight function * = 1 w(x) = (x-a)**alfa*(b-x)**beta * = 2 w(x) = (x-a)**alfa*(b-x)**beta*log(x-a) * = 3 w(x) = (x-a)**alfa*(b-x)**beta*log(b-x) * = 4 w(x) = (x-a)**alfa*(b-x)**beta*log(x-a)*log(b-x)
integer, intent(out) :: Nev

number of integrand evaluations

Calls

proc~~dqc25s~~CallsGraph proc~dqc25s quadpack_generic::dqc25s proc~dqcheb quadpack_generic::dqcheb proc~dqc25s->proc~dqcheb proc~dqk15w quadpack_generic::dqk15w proc~dqc25s->proc~dqk15w

Called by

proc~~dqc25s~~CalledByGraph proc~dqc25s quadpack_generic::dqc25s proc~dqawse quadpack_generic::dqawse proc~dqawse->proc~dqc25s proc~dqaws quadpack_generic::dqaws proc~dqaws->proc~dqawse

Source Code

    subroutine dqc25s(f, a, b, Bl, Br, Alfa, Beta, Ri, Rj, Rg, Rh, Result, Abserr, &
                      Resasc, Integr, Nev)
        implicit none

        procedure(func) :: f !! function subprogram defining the integrand f(x).
        real(wp), intent(in) :: a !! left end point of the original interval
        real(wp), intent(in) :: b !! right end point of the original interval, `b>a`
        real(wp), intent(in) :: Bl !! lower limit of integration, `bl>=a`
        real(wp), intent(in) :: Br !! upper limit of integration, `br<=b`
        real(wp), intent(in) :: Alfa !! parameter in the weight function
        real(wp), intent(in) :: Beta !! parameter in the weight function
        real(wp), intent(in) :: Ri(25) !! modified chebyshev moments for the application
                                       !! of the generalized clenshaw-curtis
                                       !! method (computed in subroutine [[dqmomo]])
        real(wp), intent(in) :: Rj(25) !! modified chebyshev moments for the application
                                       !! of the generalized clenshaw-curtis
                                       !! method (computed in subroutine [[dqmomo]])
        real(wp), intent(in) :: Rg(25) !! modified chebyshev moments for the application
                                       !! of the generalized clenshaw-curtis
                                       !! method (computed in subroutine [[dqmomo]])
        real(wp), intent(in) :: Rh(25) !! modified chebyshev moments for the application
                                       !! of the generalized clenshaw-curtis
                                       !! method (computed in subroutine [[dqmomo]])
        real(wp), intent(out) :: Result !! approximation to the integral
                                        !! `result` is computed by using a generalized
                                        !! clenshaw-curtis method if `b1 = a` or `br = b`.
                                        !! in all other cases the 15-point kronrod
                                        !! rule is applied, obtained by optimal addition of
                                        !! abscissae to the 7-point gauss rule.
        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) :: Resasc !! approximation to the integral of abs(f*w-i/(b-a))
        integer, intent(in) :: Integr !! which determines the weight function
                                      !! * = 1  `w(x) = (x-a)**alfa*(b-x)**beta`
                                      !! * = 2  `w(x) = (x-a)**alfa*(b-x)**beta*log(x-a)`
                                      !! * = 3  `w(x) = (x-a)**alfa*(b-x)**beta*log(b-x)`
                                      !! * = 4  `w(x) = (x-a)**alfa*(b-x)**beta*log(x-a)*log(b-x)`
        integer, intent(out) :: Nev !! number of integrand evaluations

        real(wp) :: cheb12(13) !! coefficients of the chebyshev series expansion
                               !! of degree 12, for the function `f`, in the
                               !! interval `(bl,br)`
        real(wp) :: cheb24(25) !! coefficients of the chebyshev series expansion
                               !! of degree 24, for the function `f`, in the
                               !! interval `(bl,br)`
        real(wp) :: fval(25) !! value of the function f at the points
                             !! `(br-bl)*0.5*cos(k*pi/24)+(br+bl)*0.5`
                             !! `k = 0, ..., 24`
        real(wp) :: res12 !! approximation to the integral obtained from `cheb12`
        real(wp) :: res24 !! approximation to the integral obtained from `cheb24`
        real(wp) :: hlgth !! half-length of the interval `(bl,br)`
        real(wp) :: centr !! mid point of the interval `(bl,br)`
        integer :: k !! counter for `x`
        real(wp) :: dc, factor, fix, resabs, u
        integer :: i, isym

        real(wp), dimension(11), parameter :: x = [(cos(k*pi/24.0_wp), k=1, 11)]
            !! the vector x contains the values `cos(k*pi/24)`,
            !! `k = 1, ..., 11`, to be used for the chebyshev series
            !! expansion of `f`

        Nev = 25
        if (Bl == a .and. (Alfa /= 0.0_wp .or. Integr == 2 .or. Integr == 4)) &
            then

            ! this part of the program is executed only if a = bl.

            ! compute the chebyshev series expansion of the
            ! following function
            ! f1 = (0.5*(b+b-br-a)-0.5*(br-a)*x)**beta
            !      *f(0.5*(br-a)*x+0.5*(br+a))

            hlgth = 0.5_wp*(Br - Bl)
            centr = 0.5_wp*(Br + Bl)
            fix = b - centr
            fval(1) = 0.5_wp*f(hlgth + centr)*(fix - hlgth)**Beta
            fval(13) = f(centr)*(fix**Beta)
            fval(25) = 0.5_wp*f(centr - hlgth)*(fix + hlgth)**Beta
            do i = 2, 12
                u = hlgth*x(i - 1)
                isym = 26 - i
                fval(i) = f(u + centr)*(fix - u)**Beta
                fval(isym) = f(centr - u)*(fix + u)**Beta
            end do
            factor = hlgth**(Alfa + 1.0_wp)
            Result = 0.0_wp
            Abserr = 0.0_wp
            res12 = 0.0_wp
            res24 = 0.0_wp
            if (Integr > 2) then

                ! compute the chebyshev series expansion of the
                ! following function
                ! f4 = f1*log(0.5*(b+b-br-a)-0.5*(br-a)*x)

                fval(1) = fval(1)*log(fix - hlgth)
                fval(13) = fval(13)*log(fix)
                fval(25) = fval(25)*log(fix + hlgth)
                do i = 2, 12
                    u = hlgth*x(i - 1)
                    isym = 26 - i
                    fval(i) = fval(i)*log(fix - u)
                    fval(isym) = fval(isym)*log(fix + u)
                end do
                call dqcheb(x, fval, cheb12, cheb24)

                ! integr = 3  (or 4)

                do i = 1, 13
                    res12 = res12 + cheb12(i)*Ri(i)
                    res24 = res24 + cheb24(i)*Ri(i)
                end do
                do i = 14, 25
                    res24 = res24 + cheb24(i)*Ri(i)
                end do
                if (Integr /= 3) then

                    ! integr = 4

                    dc = log(Br - Bl)
                    Result = res24*dc
                    Abserr = abs((res24 - res12)*dc)
                    res12 = 0.0_wp
                    res24 = 0.0_wp
                    do i = 1, 13
                        res12 = res12 + cheb12(i)*Rg(i)
                        res24 = res24 + cheb24(i)*Rg(i)
                    end do
                    do i = 14, 25
                        res24 = res24 + cheb24(i)*Rg(i)
                    end do
                end if
            else
                call dqcheb(x, fval, cheb12, cheb24)

                ! integr = 1  (or 2)

                do i = 1, 13
                    res12 = res12 + cheb12(i)*Ri(i)
                    res24 = res24 + cheb24(i)*Ri(i)
                end do
                do i = 14, 25
                    res24 = res24 + cheb24(i)*Ri(i)
                end do
                if (Integr /= 1) then

                    ! integr = 2

                    dc = log(Br - Bl)
                    Result = res24*dc
                    Abserr = abs((res24 - res12)*dc)
                    res12 = 0.0_wp
                    res24 = 0.0_wp
                    do i = 1, 13
                        res12 = res12 + cheb12(i)*Rg(i)
                        res24 = res12 + cheb24(i)*Rg(i)
                    end do
                    do i = 14, 25
                        res24 = res24 + cheb24(i)*Rg(i)
                    end do
                end if
            end if
            Result = (Result + res24)*factor
            Abserr = (Abserr + abs(res24 - res12))*factor
        elseif (Br == b .and. (Beta /= 0.0_wp .or. Integr == 3 .or. Integr == 4)) then

            ! this part of the program is executed only if b = br.

            ! compute the chebyshev series expansion of the
            ! following function
            ! f2 = (0.5*(b+bl-a-a)+0.5*(b-bl)*x)**alfa
            !      *f(0.5*(b-bl)*x+0.5*(b+bl))

            hlgth = 0.5_wp*(Br - Bl)
            centr = 0.5_wp*(Br + Bl)
            fix = centr - a
            fval(1) = 0.5_wp*f(hlgth + centr)*(fix + hlgth)**Alfa
            fval(13) = f(centr)*(fix**Alfa)
            fval(25) = 0.5_wp*f(centr - hlgth)*(fix - hlgth)**Alfa
            do i = 2, 12
                u = hlgth*x(i - 1)
                isym = 26 - i
                fval(i) = f(u + centr)*(fix + u)**Alfa
                fval(isym) = f(centr - u)*(fix - u)**Alfa
            end do
            factor = hlgth**(Beta + 1.0_wp)
            Result = 0.0_wp
            Abserr = 0.0_wp
            res12 = 0.0_wp
            res24 = 0.0_wp
            if (Integr == 2 .or. Integr == 4) then

                ! compute the chebyshev series expansion of the
                ! following function
                ! f3 = f2*log(0.5*(b-bl)*x+0.5*(b+bl-a-a))

                fval(1) = fval(1)*log(hlgth + fix)
                fval(13) = fval(13)*log(fix)
                fval(25) = fval(25)*log(fix - hlgth)
                do i = 2, 12
                    u = hlgth*x(i - 1)
                    isym = 26 - i
                    fval(i) = fval(i)*log(u + fix)
                    fval(isym) = fval(isym)*log(fix - u)
                end do
                call dqcheb(x, fval, cheb12, cheb24)

                ! integr = 2  (or 4)

                do i = 1, 13
                    res12 = res12 + cheb12(i)*Rj(i)
                    res24 = res24 + cheb24(i)*Rj(i)
                end do
                do i = 14, 25
                    res24 = res24 + cheb24(i)*Rj(i)
                end do
                if (Integr /= 2) then
                    dc = log(Br - Bl)
                    Result = res24*dc
                    Abserr = abs((res24 - res12)*dc)
                    res12 = 0.0_wp
                    res24 = 0.0_wp

                    ! integr = 4

                    do i = 1, 13
                        res12 = res12 + cheb12(i)*Rh(i)
                        res24 = res24 + cheb24(i)*Rh(i)
                    end do
                    do i = 14, 25
                        res24 = res24 + cheb24(i)*Rh(i)
                    end do
                end if
            else

                ! integr = 1  (or 3)

                call dqcheb(x, fval, cheb12, cheb24)
                do i = 1, 13
                    res12 = res12 + cheb12(i)*Rj(i)
                    res24 = res24 + cheb24(i)*Rj(i)
                end do
                do i = 14, 25
                    res24 = res24 + cheb24(i)*Rj(i)
                end do
                if (Integr /= 1) then

                    ! integr = 3

                    dc = log(Br - Bl)
                    Result = res24*dc
                    Abserr = abs((res24 - res12)*dc)
                    res12 = 0.0_wp
                    res24 = 0.0_wp
                    do i = 1, 13
                        res12 = res12 + cheb12(i)*Rh(i)
                        res24 = res24 + cheb24(i)*Rh(i)
                    end do
                    do i = 14, 25
                        res24 = res24 + cheb24(i)*Rh(i)
                    end do
                end if
            end if
            Result = (Result + res24)*factor
            Abserr = (Abserr + abs(res24 - res12))*factor
        else

            ! if a>bl and b<br, apply the 15-point gauss-kronrod
            ! scheme.

            ! dqwgts - external function subprogram defining
            ! the four possible weight functions

            call dqk15w(f, dqwgts, a, b, Alfa, Beta, Integr, Bl, Br, Result, Abserr, &
                        resabs, Resasc)
            Nev = 15
        end if

    end subroutine dqc25s