Artifact f4d620318e9aba9ccbd6e7973dfc8dbb1815c139:
# classic_polyns.tcl --
# Implement procedures for the classic orthogonal polynomials
#
package require math::polynomials
namespace eval ::math::special {
if {[info commands addPolyn] == {} } {
namespace import ::math::polynomials::*
}
}
# legendre --
# Return the nth degree Legendre polynomial
#
# Arguments:
# n The degree of the polynomial
# Result:
# Polynomial definition
#
proc ::math::special::legendre {n} {
if { ! [string is integer -strict $n] || $n < 0 } {
return -code error "Degree must be a non-negative integer"
}
set pnm1 [polynomial 1.0]
set pn [polynomial {0.0 1.0}]
if { $n == 0 } {return $pnm1}
if { $n == 1 } {return $pn}
set degree 1
while { $degree < $n } {
set an [expr {(2.0*$degree+1.0)/($degree+1.0)}]
set bn 0.0
set cn [expr {$degree/($degree+1.0)}]
set factor_n [polynomial [list $bn $an]]
set term_nm1 [multPolyn $pnm1 [expr {-1.0*$cn}]]
set term_n [multPolyn $factor_n $pn]
set pnp1 [addPolyn $term_n $term_nm1]
set pnm1 $pn
set pn $pnp1
incr degree
}
return $pnp1
}
# chebyshev --
# Return the nth degree Chebeyshev polynomial of the first kind
#
# Arguments:
# n The degree of the polynomial
# Result:
# Polynomial definition
#
proc ::math::special::chebyshev {n} {
if { ! [string is integer -strict $n] || $n < 0 } {
return -code error "Degree must be a non-negative integer"
}
set pnm1 [polynomial 1.0]
set pn [polynomial {0.0 1.0}]
if { $n == 0 } {return $pnm1}
if { $n == 1 } {return $pn}
set degree 1
while { $degree < $n } {
set an 2.0
set bn 0.0
set cn 1.0
set factor_n [polynomial [list $bn $an]]
set term_nm1 [multPolyn $pnm1 [expr {-1.0*$cn}]]
set term_n [multPolyn $factor_n $pn]
set pnp1 [addPolyn $term_n $term_nm1]
set pnm1 $pn
set pn $pnp1
incr degree
}
return $pnp1
}
# laguerre --
# Return the nth degree Laguerre polynomial with parameter alpha
#
# Arguments:
# alpha The parameter for the polynomial
# n The degree of the polynomial
# Result:
# Polynomial definition
#
proc ::math::special::laguerre {alpha n} {
if { ! [string is double -strict $alpha] } {
return -code error "Parameter must be a double"
}
if { ! [string is integer -strict $n] || $n < 0 } {
return -code error "Degree must be a non-negative integer"
}
set pnm1 [polynomial 1.0]
set pn [polynomial [list [expr {1.0-$alpha}] -1.0]]
if { $n == 0 } {return $pnm1}
if { $n == 1 } {return $pn}
set degree 1
while { $degree < $n } {
set an [expr {-1.0/($degree+1.0)}]
set bn [expr {(2.0*$degree+$alpha+1)/($degree+1.0)}]
set cn [expr {($degree+$alpha)/($degree+1.0)}]
set factor_n [polynomial [list $bn $an]]
set term_nm1 [multPolyn $pnm1 [expr {-1.0*$cn}]]
set term_n [multPolyn $factor_n $pn]
set pnp1 [addPolyn $term_n $term_nm1]
set pnm1 $pn
set pn $pnp1
incr degree
}
return $pnp1
}
# hermite --
# Return the nth degree Hermite polynomial
#
# Arguments:
# n The degree of the polynomial
# Result:
# Polynomial definition
#
proc ::math::special::hermite {n} {
if { ! [string is integer -strict $n] || $n < 0 } {
return -code error "Degree must be a non-negative integer"
}
set pnm1 [polynomial 1.0]
set pn [polynomial {0.0 2.0}]
if { $n == 0 } {return $pnm1}
if { $n == 1 } {return $pn}
set degree 1
while { $degree < $n } {
set an 2.0
set bn 0.0
set cn [expr {2.0*$degree}]
set factor_n [polynomial [list $bn $an]]
set term_n [multPolyn $factor_n $pn]
set term_nm1 [multPolyn $pnm1 [expr {-1.0*$cn}]]
set pnp1 [addPolyn $term_n $term_nm1]
set pnm1 $pn
set pn $pnp1
incr degree
}
return $pnp1
}
# some tests --
#
if { 0 } {
set tcl_precision 17
puts "Legendre:"
foreach n {0 1 2 3 4} {
puts [::math::special::legendre $n]
}
puts "Chebyshev:"
foreach n {0 1 2 3 4} {
puts [::math::special::chebyshev $n]
}
puts "Laguerre (alpha=0):"
foreach n {0 1 2 3 4} {
puts [::math::special::laguerre 0.0 $n]
}
puts "Laguerre (alpha=1):"
foreach n {0 1 2 3 4} {
puts [::math::special::laguerre 1.0 $n]
}
puts "Hermite:"
foreach n {0 1 2 3 4} {
puts [::math::special::hermite $n]
}
}