dqmomo Subroutine

public subroutine dqmomo(Alfa, Beta, Ri, Rj, Rg, Rh, Integr)

1D integration of k-th degree Chebyshev polynomial times a function with singularities

this routine computes modified chebsyshev moments. the k-th modified chebyshev moment is defined as the integral over (-1,1) of w(x)*t(k,x), where t(k,x) is the chebyshev polynomial of degree k.

History

  • QUADPACK: date written 820101, revision date 830518 (yymmdd).

Arguments

Type IntentOptional Attributes Name
real(kind=wp), intent(in) :: Alfa

parameter in the weight function w(x), alfa>(-1)
real(kind=wp), intent(in) :: Beta

parameter in the weight function w(x), beta>(-1)
real(kind=wp), intent(out) :: Ri(25)

i(k) is the integral over (-1,1) of (1+x)**alfa*t(k-1,x), k = 1, ..., 25.
real(kind=wp), intent(out) :: Rj(25)

rj(k) is the integral over (-1,1) of (1-x)**beta*t(k-1,x), k = 1, ..., 25.
real(kind=wp), intent(out) :: Rg(25)

rg(k) is the integral over (-1,1) of (1+x)**alfa*log((1+x)/2)*t(k-1,x), k = 1, ..., 25.
real(kind=wp), intent(out) :: Rh(25)

rh(k) is the integral over (-1,1) of (1-x)**beta*log((1-x)/2)*t(k-1,x), k = 1, ..., 25.
integer, intent(in) :: Integr

input parameter indicating the modified moments to be computed:

  • integr = 1 compute ri, rj
  • integr = 2 compute ri, rj, rg
  • integr = 3 compute ri, rj, rh
  • integr = 4 compute ri, rj, rg, rh

Called by

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

Source Code

    subroutine dqmomo(Alfa, Beta, Ri, Rj, Rg, Rh, Integr)
        implicit none

        real(wp), intent(in) :: Alfa !! parameter in the weight function `w(x)`, `alfa>(-1)`
        real(wp), intent(in) :: Beta !! parameter in the weight function `w(x)`, `beta>(-1)`
        real(wp), intent(out) :: Ri(25) !! `i(k)` is the integral over (-1,1) of
                                        !! `(1+x)**alfa*t(k-1,x), k = 1, ..., 25`.
        real(wp), intent(out) :: Rj(25) !! `rj(k)` is the integral over (-1,1) of
                                        !! `(1-x)**beta*t(k-1,x), k = 1, ..., 25`.
        real(wp), intent(out) :: Rg(25) !! `rg(k)` is the integral over (-1,1) of
                                        !! `(1+x)**alfa*log((1+x)/2)*t(k-1,x), k = 1, ..., 25`.
        real(wp), intent(out) :: Rh(25) !! `rh(k)` is the integral over (-1,1) of
                                        !! `(1-x)**beta*log((1-x)/2)*t(k-1,x), k = 1, ..., 25`.
        integer, intent(in) :: Integr !! input parameter indicating the modified
                                      !! moments to be computed:
                                      !!
                                      !! * integr = 1 compute `ri`, `rj`
                                      !! * integr = 2 compute `ri`, `rj`, `rg`
                                      !! * integr = 3 compute `ri`, `rj`, `rh`
                                      !! * integr = 4 compute `ri`, `rj`, `rg`, `rh`

        real(wp) :: alfp1, alfp2, an, anm1, betp1, betp2, ralf, rbet
        integer :: i, im1

        main : block

            alfp1 = Alfa + 1.0_wp
            betp1 = Beta + 1.0_wp
            alfp2 = Alfa + 2.0_wp
            betp2 = Beta + 2.0_wp
            ralf = 2.0_wp**alfp1
            rbet = 2.0_wp**betp1

            ! compute ri, rj using a forward recurrence relation.

            Ri(1) = ralf/alfp1
            Rj(1) = rbet/betp1
            Ri(2) = Ri(1)*Alfa/alfp2
            Rj(2) = Rj(1)*Beta/betp2
            an = 2.0_wp
            anm1 = 1.0_wp
            do i = 3, 25
                Ri(i) = -(ralf + an*(an - alfp2)*Ri(i - 1))/(anm1*(an + alfp1))
                Rj(i) = -(rbet + an*(an - betp2)*Rj(i - 1))/(anm1*(an + betp1))
                anm1 = an
                an = an + 1.0_wp
            end do
            if (Integr /= 1) then
                if (Integr /= 3) then

                    ! compute rg using a forward recurrence relation.

                    Rg(1) = -Ri(1)/alfp1
                    Rg(2) = -(ralf + ralf)/(alfp2*alfp2) - Rg(1)
                    an = 2.0_wp
                    anm1 = 1.0_wp
                    im1 = 2
                    do i = 3, 25
                        Rg(i) = -(an*(an - alfp2)*Rg(im1) - an*Ri(im1) + anm1*Ri(i)) &
                                /(anm1*(an + alfp1))
                        anm1 = an
                        an = an + 1.0_wp
                        im1 = i
                    end do
                    if (Integr == 2) exit main
                end if

                ! compute rh using a forward recurrence relation.

                Rh(1) = -Rj(1)/betp1
                Rh(2) = -(rbet + rbet)/(betp2*betp2) - Rh(1)
                an = 2.0_wp
                anm1 = 1.0_wp
                im1 = 2
                do i = 3, 25
                    Rh(i) = -(an*(an - betp2)*Rh(im1) - an*Rj(im1) + anm1*Rj(i)) &
                            /(anm1*(an + betp1))
                    anm1 = an
                    an = an + 1.0_wp
                    im1 = i
                end do
                do i = 2, 25, 2
                    Rh(i) = -Rh(i)
                end do
            end if

        end block main

        do i = 2, 25, 2
            Rj(i) = -Rj(i)
        end do

    end subroutine dqmomo