LagrangeInterpolationPolys Subroutine

public subroutine LagrangeInterpolationPolys(x, N_in, xGP, wBary, L)

Computes all Lagrange functions evaluated at position x in [-1,1] For details see paper Barycentric Lagrange Interpolation by Berrut and Trefethen (SIAM 2004) Uses function ALMOSTEQUAL Algorithm 34, Kopriva book

Arguments

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

Coordinate
integer, intent(in) :: N_in

polynomial degree
real(kind=wp), intent(in) :: xGP(0:N_in)

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

Barycentric weights
real(kind=wp), intent(out) :: L(0:N_in)

Lagrange basis functions evaluated at x

Calls

proc~~lagrangeinterpolationpolys~~CallsGraph proc~lagrangeinterpolationpolys LagrangeInterpolationPolys proc~almostequal ALMOSTEQUAL proc~lagrangeinterpolationpolys->proc~almostequal

Source Code

SUBROUTINE LagrangeInterpolationPolys(x,N_in,xGP,wBary,L)
IMPLICIT NONE
!----------------------------------------------------------------------------------------------------------------------------
! INPUT/OUTPUT VARIABLES
REAL(wp), INTENT(IN)   :: x                !! Coordinate
INTEGER ,INTENT(IN)    :: N_in             !! polynomial degree
REAL(wp),INTENT(IN)    :: xGP(0:N_in)      !! Gauss point positions for the reference interval [-1,1]
REAL(wp),INTENT(IN)    :: wBary(0:N_in)    !! Barycentric weights
REAL(wp),INTENT(OUT)   :: L(0:N_in)        !! Lagrange basis functions evaluated at x
!----------------------------------------------------------------------------------------------------------------------------
! LOCAL VARIABLES
INTEGER            :: iGP
LOGICAL            :: xEqualGP         ! is x equal to a Gauss Point
REAL(wp)               :: DummySum
!============================================================================================================================
xEqualGP=.FALSE.
DO iGP=0,N_in
  L(iGP)=0.0_wp
  IF(ALMOSTEQUAL(x,xGP(iGP))) THEN
    L(iGP)=1.0_wp
    xEqualGP=.TRUE.
  END IF
END DO

! if x is equal to a Gauss point, L=(0,....,1,....0)
IF(xEqualGP) RETURN
DummySum=0.0_wp
DO iGP=0, N_in
  L(iGP)=wBary(iGP)/(x-xGP(iGP))
  DummySum=DummySum+L(iGP)
END DO

DO iGP=0,N_in
  L(iGP)=L(iGP)/DummySum
END DO
END SUBROUTINE LagrangeInterpolationPolys