LegGaussLobNodesAndWeights Subroutine

public subroutine LegGaussLobNodesAndWeights(N_in, xGP, wGP)

Starting with initial guess by Parter Relation, a Newton method is used to find the roots of the Legendre Polynomial Lder_(N_in), which are the positions of Gauss-Lobatto points. Uses qAndLEvaluation subroutine. algorithm 25, Kopriva

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~~leggausslobnodesandweights~~CallsGraph proc~leggausslobnodesandweights LegGaussLobNodesAndWeights proc~qandlevaluation qAndLEvaluation proc~leggausslobnodesandweights->proc~qandlevaluation

Source Code

SUBROUTINE LegGaussLobNodesAndWeights(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)                      :: q,qder,L         ! \f$ q=L_{N_in+1}-L_{N_in-1} \f$ ,qder is derivative, \f$ L=L_{N_in} \f$
REAL(wp)                      :: dx               ! Newton step
REAL(wp)                      :: cont1,cont2      !temporary variable for evaluation of parter nodes positions
!==================================================================================================================================
xGP(0)=-1.0_wp
xGP(N_in)= 1.0_wp
IF(PRESENT(wGP))THEN
  wGP(0)= 2.0_wp/REAL(N_in*(N_in+1),wp)
  wGP(N_in)=wGP(0)
END IF
IF(N_in.GT.1)THEN
  cont1=PP_Pi/REAL(N_in,wp) ! pi/N_in
  cont2=3.0_wp/(REAL(8*N_in,wp)*PP_Pi) ! 3/(8*N_in*pi)
  DO iGP=1,(N_in+1)/2-1 !since points are symmetric, only left side is computed
    xGP(iGP)=-cos(cont1*(REAL(iGP,wp)+0.25_wp)-cont2/(REAL(iGP,wp)+0.25_wp)) !initial guess
    ! Newton iteration
    DO iter=0,nIter
      CALL qAndLEvaluation(N_in,xGP(iGP),q,qder,L)
      dx=-q/qder
      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 LGL point:'
      xGP(iGP)=-cos(cont1*(REAL(iGP,wp)+0.25)-cont2/(REAL(iGP,wp)+0.25_wp)) !initial guess
      ! Newton iteration
      DO iter=0,nIter
        WRITE(UNIT_stdout,*)'iter,x^i',iter,xGP(iGP)     !DEBUG
        CALL qAndLEvaluation(N_in,xGP(iGP),q,qder,L)
        dx=-q/qder
        xGP(iGP)=xGP(iGP)+dx
        IF(abs(dx).LT.Tol*abs(xGP(iGP))) EXIT
      END DO ! iter
      CALL abort(__STAMP__, &
                 'ERROR: Legendre Gauss Lobatto nodes could not be computed up to desired precision. Code stopped!')
    END IF ! (iter.GT.nIter)
    CALL qAndLEvaluation(N_in,xGP(iGP),q,qder,L)
    xGP(N_in-iGP)=-xGP(iGP)
    IF(PRESENT(wGP))THEN
      wGP(iGP)=wGP(0)/(L*L)
      wGP(N_in-iGP)=wGP(iGP)
    END IF
  END DO ! iGP
END IF !(N_in.GT.1)
IF(mod(N_in,2) .EQ. 0) THEN
  xGP(N_in/2)=0.0_wp
  CALL qAndLEvaluation(N_in,xGP(N_in/2),q,qder,L)
  IF(PRESENT(wGP))wGP(N_in/2)=wGP(0)/(L*L)
END IF ! (mod(N_in,2) .EQ. 0)
END SUBROUTINE LegGaussLobNodesAndWeights