JacobiP Function

  1. submodule_DG_implementation.f90
  2. DG_implementation
  3. JacobiP

function JacobiP(x, alpha, beta, n)

Function to evaluate Jacobi Polynomial of type (with ) at points for order .

Arguments

Type IntentOptional AttributesName
real(kind=wp), intent(in), dimension(:):: x

The points at which to evaluate the Jacobi Polynomial.

real(kind=wp), intent(in) :: alpha

The value of .
real(kind=wp), intent(in) :: beta

The value of .
integer(kind=ip), intent(in) :: n

The order of the Jacobi polynomial.

Return Value real(kind=wp), dimension(size(x))

Returns .


Called by

proc~~jacobip~~CalledByGraph proc~jacobip JacobiP proc~vandermonde1d Vandermonde1D proc~vandermonde1d->proc~jacobip proc~gradjacobip GradJacobiP proc~gradjacobip->proc~jacobip proc~gaussweigths GaussWeigths proc~gaussweigths->proc~jacobip proc~gradvandermonde1d GradVandermonde1D proc~gradvandermonde1d->proc~gradjacobip proc~init_ref_element init_ref_element proc~init_ref_element->proc~vandermonde1d proc~init_ref_element->proc~gaussweigths proc~dmatrix1d Dmatrix1D proc~init_ref_element->proc~dmatrix1d proc~dmatrix1d->proc~gradvandermonde1d interface~init_ref_element init_ref_element interface~init_ref_element->proc~init_ref_element interface~ref_element ref_element interface~ref_element->interface~init_ref_element proc~scal_schw_init scal_schw_init proc~scal_schw_init->interface~ref_element program~test test program~test->interface~ref_element interface~scal_schw_init scal_schw_init interface~scal_schw_init->proc~scal_schw_init

Contents

Source Code

JacobiP

Source Code

    function JacobiP(x, alpha, beta, n)
    !! Function to evaluate Jacobi Polynomial \(P^{(\alpha,\beta)}_n\)
    !! of type \((\alpha,\beta)>-1\) (with \(\alpha+\beta\neq -1\)) at points
    !! \(x\) for order \(n\).
      real(wp), dimension(:), intent(in) :: x
      !! The points \(x\in[-1:1]\) at which to evaluate the Jacobi Polynomial. 
      real(wp), intent(in) :: alpha
      !! The value of \(\alpha\).
      real(wp), intent(in) :: beta
      !! The value of \(\beta\).
      integer(ip), intent(in) :: n
      !! The order of the Jacobi polynomial.
      real(wp), dimension(size(x)) :: JacobiP
      !! Returns \(P^{(\alpha,\beta)}_n(x)\).

      real(wp), dimension(n+1,size(x)) :: pl
      real(wp) :: lngammaab, lngammaa, lngammab, invsqgamma0, gamma0, gamma1
      real(wp) :: fac1, fac2
      real(wp) :: aold, anew, bnew, h1, ireal, irealp1
      integer(ip) :: i

      pl = 0.0_wp

      lngammaab = log_gamma(alpha+beta+1.0_wp)
      lngammaa  = log_gamma(alpha+1.0_wp)
      lngammab  = log_gamma(beta+1.0_wp)

      invsqgamma0 = 2.0_wp**(alpha+beta+1.0_wp)/(alpha+beta+1.0_wp)* &
                             exp(lngammaa+lngammab-lngammaab)
      gamma0 = 1.0_wp/sqrt(invsqgamma0)

      pl(1,:) = gamma0

      if (n==0) then
        JacobiP(:) = pl(1,:)
        return
      end if

      gamma1 = 0.5_wp*sqrt((alpha+beta+3.0_wp) &
                           /((alpha+1.0_wp)*(beta+1.0_wp)))*gamma0
      fac1 = (alpha+beta+2.0_wp)
      fac2 = (alpha-beta)
      pl(2,:) = gamma1 * ( fac1*x(:) + fac2 )

      if (n==1) then
        JacobiP(:) = pl(2,:)
        return
      end if

      aold = 2.0_wp / (2.0_wp+alpha+beta) * &
             sqrt ( (1.0_wp+alpha)*(1.0_wp+beta) / (3.0_wp+alpha+beta) )

      do i = 1, n-1
        ireal = real(i,wp)
        irealp1 = ireal+1.0_wp
        h1 = 2.0_wp*ireal+alpha+beta
        anew = 2.0_wp/(h1+2.0_wp)* &
               sqrt( irealp1*(irealp1+alpha+beta) * &
                     (irealp1+alpha) * (irealp1+beta) / &
                     (h1+1.0_wp)/(h1+3.0_wp) )
        bnew = - (alpha**2-beta**2) / (h1*(h1+2.0_wp) )
        pl(i+2,:) = 1.0_wp / anew * ( -aold*pl(i,:) + (x(:)-bnew)*pl(i+1,:) )
        aold = anew
      end do
      JacobiP(:) = pl(n+1,:)
      return
    end function JacobiP