JacobiTransform

Discrete Jacobi transform.

Interface

INTERFACE
  MODULE PURE FUNCTION JacobiTransform(n, alpha, beta, coeff, x, w, &
    &  quadType) RESULT(ans)
    INTEGER(I4B), INTENT(IN) :: n
    !! order of jacobi polynomial
    REAL(DFP), INTENT(IN) :: alpha
    !! alpha of Jacobi polynomial > -1.0_DFP
    REAL(DFP), INTENT(IN) :: beta
    !! beta of Jacobi polynomial > -1.0_DFP
    REAL(DFP), INTENT(IN) :: coeff(0:n)
    !! nodal value (at quad points)
    REAL(DFP), INTENT(IN) :: x(0:n)
    !! quadrature points
    REAL(DFP), INTENT(IN) :: w(0:n)
    !! weights
    INTEGER(I4B), INTENT(IN) :: quadType
    !! Quadrature type, Gauss, GaussLobatto, GaussRadau, GaussRadauLeft
    !! GaussRadauRight
    REAL(DFP) :: ans(0:n)
    !! modal values  or coefficients
  END FUNCTION JacobiTransform
END INTERFACE

Interface 2

INTERFACE
  MODULE PURE FUNCTION JacobiTransform(n, alpha, beta, coeff, x, w, &
    & quadType) RESULT(ans)
    INTEGER(I4B), INTENT(IN) :: n
    !! order of polynomial
    REAL(DFP), INTENT(IN) :: alpha
    !! alpha of Jacobi polynomial > -1.0_DFP
    REAL(DFP), INTENT(IN) :: beta
    !! beta of Jacobi polynomial > -1.0_DFP
    REAL(DFP), INTENT(IN) :: coeff(0:, 1:)
    !! nodal value (at quad points)
    REAL(DFP), INTENT(IN) :: x(0:n)
    !! quadrature points
    REAL(DFP), INTENT(IN) :: w(0:n)
    !! weights
    INTEGER(I4B), INTENT(IN) :: quadType
    !! Quadrature type, Gauss, GaussLobatto, GaussRadau, GaussRadauLeft
    !! GaussRadauRight
    REAL(DFP) :: ans(0:n, 1:SIZE(coeff, 2))
    !! modal values  or coefficients for each column of val
  END FUNCTION JacobiTransform
END INTERFACE

Interface 3

INTERFACE
  MODULE FUNCTION JacobiTransform(n, alpha, beta, f, quadType) &
    & RESULT(ans)
    INTEGER(I4B), INTENT(IN) :: n
    !! order of jacobi polynomial
    REAL(DFP), INTENT(IN) :: alpha
    !! alpha of Jacobi polynomial > -1.0_DFP
    REAL(DFP), INTENT(IN) :: beta
    !! beta of Jacobi polynomial > -1.0_DFP
    PROCEDURE(iface_1DFunction), POINTER, INTENT(IN) :: f
    !! 1D space function
    INTEGER(I4B), INTENT(IN) :: quadType
    !! Quadrature type, Gauss, GaussLobatto, GaussRadau, GaussRadauLeft
    !! GaussRadauRight
    REAL(DFP) :: ans(0:n)
    !! modal values  or coefficients
  END FUNCTION JacobiTransform
END INTERFACE

Examples

️܀ See example

This example shows the usage of JacobiTransform method which is defined in [[JacobiPolynomialUtility]] MODULE. This routine performs the forward Jacobi transform.

In this example $\alpha=\beta=0.0$ (that is, proportional to Legendre polynomial)

program main
  use easifembase
  use easifemclasses
  implicit none
  integer( i4b ) :: n
  real(dfp), allocatable :: uhat(:), u( : ), pt( : ), wt(:)
  real( dfp ), parameter :: alpha=0.0_DFP, beta=0.0_DFP, tol=1.0E-10
  type(string) :: astr
  integer( i4b ), parameter :: quadType = GaussLobatto
  type(CSVFile_) :: afile
  character( len=* ), parameter :: fname = "./results/test15"
  type(PLPlot_) :: aplot

  n = 10
  call reallocate( pt, n+1, wt, n+1 )
  call JacobiQuadrature( n=n+1, alpha=alpha, beta=beta, &
    & pt=pt, wt=wt, quadType=quadType )
  u = JacobiEval(n=5_I4B, alpha=alpha, beta=beta, &
    & x=pt)
  uhat = JacobiTransform(n=n, alpha=alpha, &
    & beta=beta, coeff=u, &
    & x=pt, w=wt, quadType=quadType)
  uhat( 6 ) = uhat( 6 ) - 1.0_dFP
  CALL ok( SOFTEQ( NORM2(uhat), 0.0_DFP, tol), "test"  )

  n = 10
  call reallocate( pt, n+1, wt, n+1 )
  call JacobiQuadrature( n=n+1, alpha=alpha, beta=beta, &
    & pt=pt, wt=wt, quadType=quadType )
  u = SIN(4.0_DFP * pi * pt)
  uhat = JacobiTransform(n=n, alpha=alpha, &
    & beta=beta, coeff=u, x=pt, w=wt, quadType=quadType)

  CALL afile%initiate(filename=fname // ".csv", &
    & status="NEW", action="WRITE")
  CALL afile%open()
  CALL afile%write(val="n="//tostring(n))
  CALL afile%write(val=uhat)
  CALL afile%write()

  CALL aplot%Initiate()
  CALL aplot%Set( &
    & device="svg", &
    & filename=fname//"-%n.svg")
  CALL aplot%figure()
  CALL aplot%subplot(1,1,1)
  pt = arange(0,n,1)
  CALL aplot%setXYLim([-10.0_DFP, 50.0_DFP], [ -3.0_DFP, 3.0_DFP])
  CALL aplot%setTicks()

  CALL aplot%plot2D(x=pt,y=uhat, &
    & pointType=PS_DOT, &
    & pointSize=4.0_DFP, &
    & pointColor="b", &
    & lineWidth=0.0_DFP)
  !CALL aplot%setLabels("","","")
  !CALL aplot%Show()

  n = 20
  call reallocate( pt, n+1, wt, n+1 )
  call JacobiQuadrature( n=n+1, alpha=alpha, beta=beta, &
    & pt=pt, wt=wt, quadType=quadType )
  u = SIN(4.0_DFP * pi * pt)
  uhat = JacobiTransform(n=n, alpha=alpha, beta=beta,&
    & coeff=u, x=pt, w=wt, quadType=quadType)
  CALL afile%write(val="n="//tostring(n))
  CALL afile%write(val=uhat)
  CALL afile%write()

  pt = arange(0,n,1)
  CALL aplot%plot2D(x=pt,y=uhat, &
    & pointType=PS_DOT, &
    & pointSize=4.0_DFP, &
    & pointColor="r", &
    & lineWidth=0.0_DFP)
  !CALL aplot%setLabels("N","uhat","")
  !CALL aplot%Show()

  n = 30
  call reallocate( pt, n+1, wt, n+1 )
  call JacobiQuadrature( n=n+1, alpha=alpha, beta=beta, &
    & pt=pt, wt=wt, quadType=quadType )
  u = SIN(4.0_DFP * pi * pt)
  uhat = JacobiTransform(n=n, alpha=alpha, beta=beta, &
    & coeff=u, x=pt, w=wt, quadType=quadType)
  CALL afile%write(val="n="//tostring(n))
  CALL afile%write(val=uhat)
  CALL afile%write()

  pt = arange(0,n,1)
  CALL aplot%plot2D(x=pt,y=uhat, &
    & pointType=PS_DOT, &
    & pointSize=10.0_DFP, &
    & pointColor="k", &
    & lineWidth=0.0_DFP)
  CALL aplot%setLabels("N","uhat","")
  CALL aplot%Show()

  n = 40
  call reallocate( pt, n+1, wt, n+1 )
  call JacobiQuadrature( n=n+1, alpha=alpha, beta=beta, &
    & pt=pt, wt=wt, quadType=quadType )
  u = SIN(4.0_DFP * pi * pt)
  uhat = JacobiTransform(n=n, alpha=alpha, beta=beta, &
    & coeff=u, x=pt, w=wt, quadType=quadType)
  CALL afile%write(val="n="//tostring(n))
  CALL afile%write(val=uhat)
  CALL afile%write()

  CALL afile%deallocate()

end program main

↢

️܀ See example

This example shows the usage of JacobiTransform method which is defined in [[JacobiPolynomialUtility]] MODULE. This routine performs the forward Jacobi transform of column data.

In this example $\alpha=\beta=0.0$ (that is, proportional to Legendre polynomial)

program main
  use easifembase
  implicit none
  integer( i4b ) :: n
  real(dfp), allocatable :: uhat(:, :), u( :, : ), pt( : ), wt(:)
  real( dfp ), parameter :: alpha=0.0_DFP, beta=0.0_DFP, tol=1.0E-10
  type(string) :: astr
  integer( i4b ), parameter :: quadType = GaussLobatto

  n = 10
  call reallocate( pt, n+1, wt, n+1 )
  call JacobiQuadrature( n=n+1, alpha=alpha, beta=beta, &
    & pt=pt, wt=wt, quadType=quadType )
  u = SIN(4.0_DFP * pi * pt) .COLCONCAT. SIN(4.0_DFP * pi * pt)
  uhat = JacobiTransform(n=n, alpha=alpha, beta=beta, coeff=u, &
    & x=pt, w=wt, quadType=quadType)
  CALL Display(uhat, "uhat=")

end program main

↢

️܀ See example

This example shows the usage of JacobiTransform method.

In this example $\alpha=\beta=0.0$ (that is, proportional to Legendre polynomial)

program main
  use easifembase
  use easifemclasses
  implicit none
  integer( i4b ) :: n, ii
  real(dfp), allocatable :: uhat(:), u( : ), pt( : ), wt(:), fhat(:)
  real( dfp ), parameter :: alpha=0.0_DFP, beta=0.0_DFP, tol=1.0E-10
  type(string) :: astr
  integer( i4b ), parameter :: quadType = GaussLobatto
  procedure(iface_1Dfunction), pointer :: f => NULL()

note

To make quadrature points we call JacobiQuadrature routine, and specify QuadType.

  n = 10
  call reallocate( pt, n+1, wt, n+1, u, n+1 )
  call JacobiQuadrature( n=n+1, alpha=alpha, beta=beta, &
    & pt=pt, wt=wt, quadType=quadType )

  f => func1
  !!
  do ii = 1, size(u)
    u(ii) = f(pt(ii))
  end do
  !!
  uhat = JacobiTransform(n=n, alpha=alpha, &
    & beta=beta, coeff=u, x=pt, w=wt, quadType=quadType)
  !!
  fhat = JacobiTransform(n=n, alpha=alpha, beta=beta, &
    & f=f, quadType=quadType)
  !!
  call OK(ALL(SOFTEQ(fhat, uhat, tol)), "n=10")
  !!

  contains
  pure function func1(x) result(ans)
    real(dfp), intent(in) :: x
    real(dfp) :: ans
    ans = SIN(4.0_DFP * pi * x)
  end function func1

end program main

↢