LegendreGaussNodesAndWeights Subroutine

public subroutine LegendreGaussNodesAndWeights(N_in, xGP, wGP)

Starting with Chebychev point positions, a Newton method is used to find the roots of the Legendre Polynomial L_(N_in+1), which are the positions of Gausspoints uses LegendrePolynomialAndDerivative subroutine

Arguments

Type IntentOptional Attributes Name
integer, intent(in) :: N_in

polynomial degree, (N_in+1) Gausspoints
real(kind=wp), intent(out) :: xGP(0:N_in)

Gauss point positions for the reference interval [-1,1]
real(kind=wp), intent(out), optional :: wGP(0:N_in)

Gauss point weights

Calls

proc~~legendregaussnodesandweights~~CallsGraph proc~legendregaussnodesandweights LegendreGaussNodesAndWeights interface~legendrepolynomialandderivative LegendrePolynomialAndDerivative proc~legendregaussnodesandweights->interface~legendrepolynomialandderivative interface~legendrepolynomialandderivative->interface~legendrepolynomialandderivative

Source Code

SUBROUTINE LegendreGaussNodesAndWeights(N_in,xGP,wGP)
!MODULES
IMPLICIT NONE
!----------------------------------------------------------------------------------------------------------------------------------
! INPUT/OUTPUT VARIABLES
INTEGER,INTENT(IN)            :: N_in              !! polynomial degree, (N_in+1) Gausspoints
REAL(wp),INTENT(OUT)          :: xGP(0:N_in)       !! Gauss point positions for the reference interval [-1,1]
REAL(wp),INTENT(OUT),OPTIONAL :: wGP(0:N_in)       !! Gauss point weights
!----------------------------------------------------------------------------------------------------------------------------------
! LOCAL VARIABLES
INTEGER                   :: nIter = 10        ! max. number of newton iterations
REAL(wp)                  :: Tol   = 1.E-15_wp    ! tolerance for Newton iteration: TODO: use variable tolerance here!
INTEGER                   :: iGP,iter
REAL(wp)                  :: L_Np1,Lder_Np1    ! L_{N_in+1},Lder_{N_in+1}
REAL(wp)                  :: dx                ! Newton step
REAL(wp)                  :: cheb_tmp          ! temporary variable for evaluation of chebychev node positions
!==================================================================================================================================
IF(N_in .EQ. 0) THEN
  xGP=0.0_wp
  IF(PRESENT(wGP))wGP=2.0_wp
  RETURN
ELSEIF(N_in.EQ.1)THEN
  xGP(0)=-sqrt(1.0_wp/3.0_wp)
  xGP(N_in)=-xGP(0)
  IF(PRESENT(wGP))wGP=1.0_wp
  RETURN
ELSE ! N_in>1
  cheb_tmp=2.0_wp*ATAN(1.0_wp)/REAL(N_in+1,wp) ! pi/(2N+2)
  DO iGP=0,(N_in+1)/2-1 !since points are symmetric, only left side is computed
    xGP(iGP)=-COS(cheb_tmp*REAL(2*iGP+1,wp)) !initial guess
    ! Newton iteration
    DO iter=0,nIter
      CALL LegendrePolynomialAndDerivative(N_in+1,xGP(iGP),L_Np1,Lder_Np1)
      dx=-L_Np1/Lder_Np1
      xGP(iGP)=xGP(iGP)+dx
      IF(abs(dx).LT.Tol*abs(xGP(iGP))) EXIT
    END DO ! iter
    IF(iter.GT.nIter) THEN
      WRITE(UNIT_stdout,*) 'maximum iteration steps >10 in Newton iteration for Legendre Gausspoint'
      xGP(iGP)=-cos(cheb_tmp*REAL(2*iGP+1)) !initial guess
      ! Newton iteration
      DO iter=0,nIter
        WRITE(UNIT_stdout,*)iter,xGP(iGP)    !DEBUG
        CALL LegendrePolynomialAndDerivative(N_in+1,xGP(iGP),L_Np1,Lder_Np1)
        dx=-L_Np1/Lder_Np1
        xGP(iGP)=xGP(iGP)+dx
        IF(abs(dx).LT.Tol*abs(xGP(iGP))) EXIT
      END DO !iter
      CALL abort(__STAMP__, &
                 'ERROR: Legendre Gauss nodes could not be computed up to desired precision. Code stopped!')
    END IF ! (iter.GT.nIter)
    CALL LegendrePolynomialAndDerivative(N_in+1,xGP(iGP),L_Np1,Lder_Np1)
    xGP(N_in-iGP)=-xGP(iGP)
    IF(PRESENT(wGP))THEN
      !wGP(iGP)=2./((1.-xGP(iGP)*xGP(iGP))*Lder_Np1*Lder_Np1) !if Legendre not normalized
      wGP(iGP)=REAL(2*N_in+3,wp)/((1.0_wp-xGP(iGP)*xGP(iGP))*Lder_Np1*Lder_Np1)
      wGP(N_in-iGP)=wGP(iGP)
    END IF
  END DO !iGP
END IF ! N_in
IF(mod(N_in,2) .EQ. 0) THEN
  xGP(N_in/2)=0.0_wp
  CALL LegendrePolynomialAndDerivative(N_in+1,xGP(N_in/2),L_Np1,Lder_Np1)
  !IF(PRESENT(wGP))wGP(N_in/2)=2./(Lder_Np1*Lder_Np1) !if Legendre not normalized
  IF(PRESENT(wGP))wGP(N_in/2)=(REAL(2*N_in+3,wp))/(Lder_Np1*Lder_Np1)
END IF ! (mod(N_in,2) .EQ. 0)
END SUBROUTINE LegendreGaussNodesAndWeights