JacobiQuadrature

This routine returns the Quadrature point of Jacobi polynomial.

Here n is the number of quadrature points. Please note it is not the order of jacobi polynomial. The order is decided internally depending upon the quadType

pt and wt should be allocated outside, and length should be n.

Interface

INTERFACE
  MODULE SUBROUTINE JacobiQuadrature(n, alpha, beta, pt, wt, quadType)
    INTEGER(I4B), INTENT(IN) :: n
    !! number of quadrature points, the order will be computed as follows
    !! for quadType = Gauss, n is same as order of Jacobi polynomial
    !! for quadType = GaussRadauLeft or GaussRadauRight n is order+1
    !! for quadType = GaussLobatto, n = order+2
    REAL(DFP), INTENT(IN) :: alpha
    !! alpha of Jacobi polynomial
    REAL(DFP), INTENT(IN) :: beta
    !! beta of Jacobi polynomial
    REAL(DFP), INTENT(OUT) :: pt(n)
    !! n+1 quadrature points from 1 to n+1
    REAL(DFP), INTENT(OUT) :: wt(n)
    !! n+1 weights from 1 to n+1
    INTEGER(I4B), INTENT(IN) :: quadType
    !! Gauss
    !! GaussRadauLeft
    !! GaussRadauRight
    !! GaussLobatto
  END SUBROUTINE JacobiQuadrature
END INTERFACE

Examples

️܀ See example

This example shows the usage of JacobiQuadrature method which is defined in [[JacobiPolynomialUtility]] MODULE. This routine returns the quadrature points for Jacobi weights.

In this example $\alpha=\beta=0.0$ (that is, Legendre polynomial) By using this subroutine we can get Jacobi-Gauss, Jacobi-Gauss-Radau, Jacobi-Gauss-Lobatto quadrature points

program main
  use easifembase
  implicit none
  integer( i4b ) :: n, quadType
  real( dfp ), allocatable :: pt( : ), wt( : )
  real( dfp ), parameter :: alpha=0.0_DFP, beta=0.0_DFP
  type(string) :: msg, astr

n = 2; quadType=Gauss; call callme

See results

pt	wt
-0.57735	1
0.57735	1

n = 3; quadType=GaussRadauLeft; call callme

results

| pt      | wt      |
|---------|---------|
| -1      | 0.22222 |
| -0.2899 | 1.025   |
| 0.6899  | 0.75281 |

n = 3; quadType=GaussRadauRight; call callme

See results

pt	wt
-0.6899	0.75281
0.2899	1.025
1	0.22222

n = 4; quadType=GaussLobatto; call callme

See results

pt	wt
-1	0.16667
-0.44721	0.83333
0.44721	0.83333
1	0.16667

contains
subroutine callme
  call reallocate( pt, n, wt, n )
  call JacobiQuadrature( n=n, alpha=alpha, beta=beta, pt=pt, wt=wt, &
    & quadType=quadType )
  msg = "| pt | wt |"
  call display(msg%chars())
  astr = MdEncode( pt .COLCONCAT. wt )
  call display( astr%chars(), "" )
end subroutine callme

end program main

↢