! D01TDF Example Program Text
! Mark 28.3 Release. NAG Copyright 2022.
Module d01tdfe_mod
! D01TDF Example Program Module:
! Parameters and User-defined Routines
! .. Use Statements ..
Use nag_library, Only: nag_wp
! .. Implicit None Statement ..
Implicit None
! .. Accessibility Statements ..
Private
Public :: basic_types
! .. Parameters ..
Integer, Parameter, Public :: chebyshev1 = 2, chebyshev2 = 3, &
hermite = 6, jacobi = 4, &
laguerre = 5, legendre = 1
Contains
Subroutine basic_types(rulekind,alpha,beta,n,a,b,c,muzero)
! This procedure supplies the coefficients of the three term
! recurrence relationship for various classical orthogonal
! polynomials.
! .. Use Statements ..
Use nag_library, Only: s14aaf, x01aaf
! .. Scalar Arguments ..
Real (Kind=nag_wp), Intent (In) :: alpha, beta
Real (Kind=nag_wp), Intent (Out) :: muzero
Integer, Intent (In) :: n, rulekind
! .. Array Arguments ..
Real (Kind=nag_wp), Intent (Out) :: a(1:n), b(1:n), c(1:n)
! .. Local Scalars ..
Real (Kind=nag_wp) :: abl, mypi
Integer :: i, ifail
! .. Intrinsic Procedures ..
Intrinsic :: real, sqrt
! .. Executable Statements ..
mypi = x01aaf(mypi)
Select Case (rulekind)
Case (legendre)
muzero = 2.0_nag_wp
! P(x) in [-1, 1), w(x) = 1.0
Do i = 1, n
a(i) = (real(2*i-1,kind=nag_wp))/real(i,kind=nag_wp)
b(i) = 0.0_nag_wp
c(i) = real(i-1,kind=nag_wp)/real(i,kind=nag_wp)
End Do
Case (chebyshev1)
muzero = mypi
! c(i)= (i-1)/i
! T (x) in [-1,1], w(x) = (1-x**2)**(-0.5)
Do i = 1, n
a(i) = 2.0_nag_wp
b(i) = 0.0_nag_wp
c(i) = 1.0_nag_wp
End Do
If (n>0) Then
a(1) = 1.0_nag_wp
End If
Case (chebyshev2)
muzero = mypi/2.0_nag_wp
! u(x) in [-1, 11, W(x) = (1-x**2)** 0.5;
Do i = 1, n
a(i) = 2.0_nag_wp
b(i) = 0.0_nag_wp
c(i) = 1.0_nag_wp
End Do
Case (jacobi)
ifail = 0
muzero = 2**(alpha+beta+1)*s14aaf(alpha+1,ifail)* &
s14aaf(beta+1,ifail)/s14aaf(alpha+beta+2,ifail)
! P(alpha,beta)(x) in [-1, 11, w(x) = (1-x)**alpha*(l+x)**beta
! alpha> -1 and beta > -1
If (n>0) Then
a(1) = 0.5_nag_wp*(alpha+beta+2)
b(1) = 0.5_nag_wp*(alpha-beta)
c(1) = 0.0_nag_wp
End If
Do i = 2, n
abl = 2.0_nag_wp*real(i,kind=nag_wp)*(real(i,kind=nag_wp)+alpha+ &
beta)
a(i) = (2.0_nag_wp*real(i,kind=nag_wp)+alpha+beta-1.0_nag_wp)* &
(2.0_nag_wp*real(i,kind=nag_wp)+alpha+beta)/abl
abl = (2.0_nag_wp*real(i,kind=nag_wp)+alpha+beta-2.0_nag_wp)*abl
b(i) = (2.0_nag_wp*real(i,kind=nag_wp)+alpha+beta-1.0_nag_wp)* &
(alpha**2-beta**2)/abl
c(i) = 2.0_nag_wp*(real(i,kind=nag_wp)-1.0_nag_wp+alpha)* &
(real(i,kind=nag_wp)-1.0_nag_wp+beta)* &
(2.0_nag_wp*real(i,kind=nag_wp)+alpha+beta)/abl
End Do
Case (laguerre)
ifail = 0
muzero = s14aaf(alpha+1.0_nag_wp,ifail)
! L(alpha)(x) in [0, infinity), w(x) = exp(-x)*x**alpha,
! alpha > -1
Do i = 1, n
a(i) = -1.0_nag_wp/real(i,kind=nag_wp)
b(i) = (2.0_nag_wp*real(i,kind=nag_wp)-1.0_nag_wp+alpha)/ &
real(i,kind=nag_wp)
c(i) = (real(i,kind=nag_wp)-1.0_nag_wp+alpha)/real(i,kind=nag_wp)
End Do
Case (hermite)
muzero = sqrt(mypi)
! H(x) in (-infinity,+infinity), w(x) = exp(-x**2)
Do i = 1, n
a(i) = 2.0_nag_wp
b(i) = 0.0_nag_wp
c(i) = 2.0_nag_wp*real(i-1,kind=nag_wp)
End Do
End Select
End Subroutine basic_types
End Module d01tdfe_mod
Program d01tdfe
! .. Use Statements ..
Use d01tdfe_mod, Only: basic_types, chebyshev1, chebyshev2, hermite, &
jacobi, laguerre, legendre
Use nag_library, Only: d01tdf, nag_wp
! .. Implicit None Statement ..
Implicit None
! .. Parameters ..
Integer, Parameter :: n = 4, nout = 6
! .. Local Scalars ..
Real (Kind=nag_wp) :: alpha, beta, muzero
Integer :: ifail, j, rule
! .. Local Arrays ..
Real (Kind=nag_wp) :: a(1:n), b(1:n), c(1:n), t(1:n), &
w(1:n)
! .. Executable Statements ..
Write (nout,*) 'D01TDF Example Program Results '
Write (6,*)
rule = legendre
Do rule = 1, 6
ifail = -1
Select Case (rule)
Case (legendre)
Call basic_types(rule,alpha,beta,n,a,b,c,muzero)
Write (nout,*) 'Using the Gauss-Legendre Rule'
Call d01tdf(n,a,b,c,muzero,w,t,ifail)
If (ifail==0) Then
Write (nout,99998)
Write (nout,99997)(j,t(j),w(j),j=1,n)
Write (6,*)
End If
Case (chebyshev1)
Call basic_types(rule,alpha,beta,n,a,b,c,muzero)
Write (nout,*) 'Using the Chebyshev Rule 1'
Call d01tdf(n,a,b,c,muzero,w,t,ifail)
If (ifail==0) Then
Write (nout,99998)
Write (nout,99997)(j,t(j),w(j),j=1,n)
Write (6,*)
End If
Case (chebyshev2)
Call basic_types(rule,alpha,beta,n,a,b,c,muzero)
Write (nout,*) 'Using the Chebyshev Rule 2'
Call d01tdf(n,a,b,c,muzero,w,t,ifail)
If (ifail==0) Then
Write (nout,99998)
Write (nout,99997)(j,t(j),w(j),j=1,n)
Write (6,*)
End If
Case (jacobi)
alpha = 0.5_nag_wp
beta = 0.5_nag_wp
Call basic_types(rule,alpha,beta,n,a,b,c,muzero)
Write (nout,99999) alpha, beta
Call d01tdf(n,a,b,c,muzero,w,t,ifail)
If (ifail==0) Then
Write (nout,99998)
Write (nout,99997)(j,t(j),w(j),j=1,n)
Write (6,*)
End If
Case (laguerre)
Call basic_types(rule,alpha,beta,n,a,b,c,muzero)
Write (nout,*) 'Using the Laguerre Rule'
Call d01tdf(n,a,b,c,muzero,w,t,ifail)
If (ifail==0) Then
Write (nout,99998)
Write (nout,99997)(j,t(j),w(j),j=1,n)
Write (6,*)
End If
Case (hermite)
Call basic_types(rule,alpha,beta,n,a,b,c,muzero)
Write (nout,*) 'Using the Hermite Rule'
Call d01tdf(n,a,b,c,muzero,w,t,ifail)
If (ifail==0) Then
Write (nout,99998)
Write (nout,99997)(j,t(j),w(j),j=1,n)
Write (6,*)
End If
End Select
End Do
99999 Format (' Using the Jacobi Rule: alpha = ',F10.5,' beta = ',F10.5)
99998 Format (/,' j ',' Abscissa ',' Weight')
99997 Format (I8,D25.15,D25.15)
End Program d01tdfe