JacobiGaussRadauQuadrature

This routine returns the $n+1$ Quadrature points and weights.

The Gauss-Radau quadrature points consists one of the end points denoted by $a$. So $a$ can be $\pm 1$. The remaining $n$ points are internal to $(-1, +1)$, and they are n-zeros of Jacobi polynomial of order n with respect to the following weight.

• $(1-x)^{\alpha} (1+x)^{\beta} (x+1)$ if $a=-1$.
• $(1-x)^{\alpha} (1+x)^{\beta} (1-x)$ if $a=+1$.

Here n is the order of Jacobi polynomial.

If $a=1$ then n+1 quadrature point will be +1 If $a=-1$ then 1st quadrature point will be -1

INTERFACE
  MODULE SUBROUTINE JacobiGaussRadauQuadrature(a, n, alpha, beta, pt, wt)
    REAL(DFP), INTENT(IN) :: a
    !! the value of one of the end points
    !! it should be either -1 or +1
    INTEGER(I4B), INTENT(IN) :: n
    !! order of jacobi polynomial
    REAL(DFP), INTENT(IN) :: alpha
    !! alpha of Jacobi polynomial
    REAL(DFP), INTENT(IN) :: beta
    !! beta of Jacobi polynomial
    REAL(DFP), INTENT(OUT) :: pt(:)
    !! n+1 quadrature points from 1 to n+1
    REAL(DFP), OPTIONAL, INTENT(OUT) :: wt(:)
    !! n+1 weights from 1 to n+1
  END SUBROUTINE JacobiGaussRadauQuadrature
END INTERFACE

Examples

️܀ See example

This example shows the usage of JacobiGaussRadauQuadrature method. This routine returns the quadrature points for Jacobi weights.

In this example $\alpha=\beta=0.0$ (that is, Legendre polynomial) Note that this routine returns n+1 quadrature points and one of the points will be end point.

program main
  use easifembase
  implicit none
  integer( i4b ) :: n
  real( dfp ), allocatable :: pt( : ), wt( : )
  real( dfp ), parameter :: alpha=0.0, beta=0.0, left=-1.0, right=1.0
  type(string) :: msg, astr
  real( dfp ) :: a

order = 1

  n = 1; a=left; call callme

See results

Jacobi Gauss Radau points, n+1 = 2 alpha=0 beta=0

pt	wt
-1	0.5
0.33333	1.5

order=2.

  n = 2; a=left; call callme

See results

Jacobi Gauss Radau points, n+1 = 3 alpha=0 beta=0

pt	wt
-1	0.22222
-0.2899	1.025
0.6899	0.75281

order=1.

  n = 1; a=right; call callme

See results

Jacobi Gauss Radau points, n+1 = 2 alpha=0 beta=0

pt	wt
-0.33333	1.5
1	0.5

order = 2

  n = 2; a=right; call callme

See results

Jacobi Gauss Radau points, n+1 = 3 alpha=0 beta=0

pt	wt
-0.6899	0.75281
0.2899	1.025
1	0.22222

  contains
  subroutine callme
    call reallocate( pt, n+1, wt, n+1 )
    call JacobiGaussRadauQuadrature( a=a, n=n, &
      & alpha=alpha, beta=beta, pt=pt, wt=wt )
    msg="Jacobi Gauss Radau points, n+1 = " &
        & // tostring( n+1 ) // " alpha="//tostring( alpha ) // &
        & " beta="//tostring( beta )
    call display(msg%chars())
    astr = MdEncode( pt .COLCONCAT. wt )
    call display( astr%chars(), "" )
  end subroutine callme

end program main

↢