Fortran-90 driver and function for ChebyU
program
ChebyU_
!
Function Name: ChebyU
!
Chebyshev Polynomial of second kind, order n
double precision ChebyU,x,eps
integer n
write(*,*) "Chebyshev
Polynomial (U) Tests"
x = 0.0d00
do while ( abs(x) <= 1.0d00 )
write(*,*)
write(*,*) "Input n,x,eps (|x|>1 to end):"
read(*,*) n,x,eps
if ( abs(x) > 1.0d00 ) then
write(*,*) &
"End Chebyshev Polynomial (U) Tests"
stop
end if
write(*,*) ChebyU(n,x,eps)
end do
end program ChebyU_
function
ChebyU(n,x,eps)
! Chebyshev Polynomial of second kind, order n
double precision ChebyU,x,eps, &
xNearOne,nPone,nTerm,sign,theta
integer n
if ( n < 1 ) then
ChebyU = 1.0d00; return
end if
if ( n < 2 ) then
ChebyU = 2.0d00*x; return !
U1
end if
if ( n < 3 ) then
ChebyU = -1.0d00+4.0d00*x*x; return !
U2
end if
if ( n < 4 ) then
ChebyU = 4.0d00*x*(-1.0d00+2.0d00*x*x)
return !
U3
end if
! For n >= 4 use trig formulas
!
First check if |x| near 1
xNearOne = 1.0d00 - x*x
nPone = n+1
nTerm = n*(n+2)
if ( xNearOne <
&
sqrt(120.0d00*eps)/(nPone*nPone) ) then
!
Within danger zone
for accuracy eps
!
Sign negative if n
is odd & x is negative:
!
else sign is positive
if ( (n>2*(n/2)).and.(x<0.0d00) ) then
sign = - 1.0d00
else
sign = + 1.0d00
end if
ChebyU = sign*nPone* &
(1.0d00 - nTerm*xNearOne/6.0d00)
return
end if
!
Outside danger zone for fractional accuracy eps
theta = acos(x)
ChebyU = sin(nPone*theta)/sin(theta)
end function ChebyU