1# classic_polyns.tcl -- 2# Implement procedures for the classic orthogonal polynomials 3# 4package require math::polynomials 5 6namespace eval ::math::special { 7 if {[info commands addPolyn] == {} } { 8 namespace import ::math::polynomials::* 9 } 10} 11 12 13# legendre -- 14# Return the nth degree Legendre polynomial 15# 16# Arguments: 17# n The degree of the polynomial 18# Result: 19# Polynomial definition 20# 21proc ::math::special::legendre {n} { 22 if { ! [string is integer -strict $n] || $n < 0 } { 23 return -code error "Degree must be a non-negative integer" 24 } 25 26 set pnm1 [polynomial 1.0] 27 set pn [polynomial {0.0 1.0}] 28 29 if { $n == 0 } {return $pnm1} 30 if { $n == 1 } {return $pn} 31 32 set degree 1 33 while { $degree < $n } { 34 set an [expr {(2.0*$degree+1.0)/($degree+1.0)}] 35 set bn 0.0 36 set cn [expr {$degree/($degree+1.0)}] 37 set factor_n [polynomial [list $bn $an]] 38 set term_nm1 [multPolyn $pnm1 [expr {-1.0*$cn}]] 39 set term_n [multPolyn $factor_n $pn] 40 set pnp1 [addPolyn $term_n $term_nm1] 41 42 set pnm1 $pn 43 set pn $pnp1 44 incr degree 45 } 46 47 return $pnp1 48} 49 50# chebyshev -- 51# Return the nth degree Chebeyshev polynomial of the first kind 52# 53# Arguments: 54# n The degree of the polynomial 55# Result: 56# Polynomial definition 57# 58proc ::math::special::chebyshev {n} { 59 if { ! [string is integer -strict $n] || $n < 0 } { 60 return -code error "Degree must be a non-negative integer" 61 } 62 63 set pnm1 [polynomial 1.0] 64 set pn [polynomial {0.0 1.0}] 65 66 if { $n == 0 } {return $pnm1} 67 if { $n == 1 } {return $pn} 68 69 set degree 1 70 while { $degree < $n } { 71 set an 2.0 72 set bn 0.0 73 set cn 1.0 74 set factor_n [polynomial [list $bn $an]] 75 set term_nm1 [multPolyn $pnm1 [expr {-1.0*$cn}]] 76 set term_n [multPolyn $factor_n $pn] 77 set pnp1 [addPolyn $term_n $term_nm1] 78 79 set pnm1 $pn 80 set pn $pnp1 81 incr degree 82 } 83 84 return $pnp1 85} 86 87# laguerre -- 88# Return the nth degree Laguerre polynomial with parameter alpha 89# 90# Arguments: 91# alpha The parameter for the polynomial 92# n The degree of the polynomial 93# Result: 94# Polynomial definition 95# 96proc ::math::special::laguerre {alpha n} { 97 if { ! [string is double -strict $alpha] } { 98 return -code error "Parameter must be a double" 99 } 100 if { ! [string is integer -strict $n] || $n < 0 } { 101 return -code error "Degree must be a non-negative integer" 102 } 103 104 set pnm1 [polynomial 1.0] 105 set pn [polynomial [list [expr {1.0-$alpha}] -1.0]] 106 107 if { $n == 0 } {return $pnm1} 108 if { $n == 1 } {return $pn} 109 110 set degree 1 111 while { $degree < $n } { 112 set an [expr {-1.0/($degree+1.0)}] 113 set bn [expr {(2.0*$degree+$alpha+1)/($degree+1.0)}] 114 set cn [expr {($degree+$alpha)/($degree+1.0)}] 115 set factor_n [polynomial [list $bn $an]] 116 set term_nm1 [multPolyn $pnm1 [expr {-1.0*$cn}]] 117 set term_n [multPolyn $factor_n $pn] 118 set pnp1 [addPolyn $term_n $term_nm1] 119 120 set pnm1 $pn 121 set pn $pnp1 122 incr degree 123 } 124 125 return $pnp1 126} 127 128# hermite -- 129# Return the nth degree Hermite polynomial 130# 131# Arguments: 132# n The degree of the polynomial 133# Result: 134# Polynomial definition 135# 136proc ::math::special::hermite {n} { 137 if { ! [string is integer -strict $n] || $n < 0 } { 138 return -code error "Degree must be a non-negative integer" 139 } 140 141 set pnm1 [polynomial 1.0] 142 set pn [polynomial {0.0 2.0}] 143 144 if { $n == 0 } {return $pnm1} 145 if { $n == 1 } {return $pn} 146 147 set degree 1 148 while { $degree < $n } { 149 set an 2.0 150 set bn 0.0 151 set cn [expr {2.0*$degree}] 152 set factor_n [polynomial [list $bn $an]] 153 set term_n [multPolyn $factor_n $pn] 154 set term_nm1 [multPolyn $pnm1 [expr {-1.0*$cn}]] 155 set pnp1 [addPolyn $term_n $term_nm1] 156 157 set pnm1 $pn 158 set pn $pnp1 159 incr degree 160 } 161 162 return $pnp1 163} 164 165# some tests -- 166# 167if { 0 } { 168set tcl_precision 17 169 170puts "Legendre:" 171foreach n {0 1 2 3 4} { 172 puts [::math::special::legendre $n] 173} 174 175puts "Chebyshev:" 176foreach n {0 1 2 3 4} { 177 puts [::math::special::chebyshev $n] 178} 179 180puts "Laguerre (alpha=0):" 181foreach n {0 1 2 3 4} { 182 puts [::math::special::laguerre 0.0 $n] 183} 184puts "Laguerre (alpha=1):" 185foreach n {0 1 2 3 4} { 186 puts [::math::special::laguerre 1.0 $n] 187} 188 189puts "Hermite:" 190foreach n {0 1 2 3 4} { 191 puts [::math::special::hermite $n] 192} 193} 194