LegendreGauss Subroutine

  1. SELF_Quadrature.f90
  2. SELF_Quadrature
  3. LegendreGauss

private subroutine LegendreGauss(N, nodes, weights)

Arguments

TypeIntentOptionalAttributesName
integer :: N
real(kind=real64) :: nodes(0:N)
real(kind=real64) :: weights(0:N)

Calls

proc~~legendregauss~~CallsGraph proc~legendregauss LegendreGauss proc~legendrepolynomial LegendrePolynomial proc~legendregauss->proc~legendrepolynomial

Called by

proc~~legendregauss~~CalledByGraph proc~legendregauss LegendreGauss proc~legendrequadrature LegendreQuadrature proc~legendrequadrature->proc~legendregauss proc~init_lagrange_t Init_Lagrange_t proc~init_lagrange_t->proc~legendrequadrature proc~init_lagrange~2 Init_Lagrange proc~init_lagrange~2->proc~legendrequadrature proc~init_lagrange Init_Lagrange proc~init_lagrange->proc~legendrequadrature

Contents

Source Code

LegendreGauss

Source Code

  subroutine LegendreGauss(N,nodes,weights)
    implicit none
    integer    :: N
    real(real64) :: nodes(0:N)
    real(real64) :: weights(0:N)
    ! Local
    real(real64) :: nodes_local(0:N)
    real(real64) :: weights_local(0:N)
    real(real64) :: lN1,dlN1
    real(real64) :: delta
    integer  :: j,kIt

    if(N == 0) then

      nodes_local(0) = 0.0_real64
      weights_local(0) = 2.0_real64

    elseif(N == 1) then

      nodes_local(0) = -sqrt(1.0_real64/3.0_real64)
      weights_local(0) = 1.0_real64
      nodes_local(1) = -nodes(0)
      weights_local(1) = weights(0)

    else

      do j = 0,((N+1)/2)

        nodes_local(j) = -cos((2.0_real64*real(j,real64)+1.0_real64)*pi/(2.0_real64*real(N,real64)+1.0_real64))

        do kIt = 1,newtonMax

          call LegendrePolynomial(N+1,nodes_local(j),lN1,dlN1)
          delta = -lN1/dlN1
          nodes_local(j) = nodes_local(j)+delta
          if(abs(delta) <= TOL*nodes_local(j)) exit

        enddo

        call LegendrePolynomial(N+1,nodes_local(j),lN1,dlN1)
        weights_local(j) = 2.0_real64/((1.0_real64-nodes_local(j)*nodes_local(j))*dlN1*dlN1)
        weights_local(N-j) = weights_local(j)
        nodes_local(N-j) = -nodes_local(j)

      enddo

    endif

    if(mod(real(N,real64),2.0_real64) == 0.0_real64) then

      call LegendrePolynomial(N+1,0.0_real64,lN1,dlN1)
      nodes_local(N/2) = 0.0_real64
      weights_local(N/2) = 2.0/(dlN1*dlN1)

    endif

    do j = 0,N
      nodes(j) = real(nodes_local(j),real64)
      weights(j) = real(weights_local(j),real64)
    enddo

  endsubroutine LegendreGauss