AM (15 july 2004) The code below illustrates graphically how approximating functions with Chebyshev polynomials works, it also forms a nice little package of mathematical methods. Implemented by
kbk, put on the Wiki by me.
namespace eval ::math::approximate {
namespace export cheby_fit
}
#----------------------------------------------------------------------
#
# cheby_eval --
#
# Approximates a function with a truncated Chebyshev series
#
# Parameters:
# coef - List of coefficients of the Chebyshev polynomials
# a, b - Boundaries of the interval over which the Chebyshev
# fit was computed.
# x - Abscissa at which to evaluate the function
#
# Results:
# Returns an approximation to the function evaluated at x.
#
# Side effects:
# None.
#
#----------------------------------------------------------------------
proc ::math::approximate::cheby_eval { coef a b x } {
if { ( $x - $a ) * ( $x - $b ) > 1e-8 } {
return -code error "domain error: x not in interval"
}
set u [expr { ( $x + $x - $a - $b ) / ( $b - $a ) }]
set twou [expr { $u + $u }]
set v2 0.
set v1 0.
set j [llength $coef]
while { $j >= 2 } {
incr j -1
foreach { v2 v1 } \
[list $v1 [expr { $twou * $v1 - $v2 + [lindex $coef $j] }]] \
break
}
return [expr { $u * $v1 - $v2 + [lindex $coef 0] }]
}
#----------------------------------------------------------------------
#
# cheby_fit --
#
# Fits a function with a truncated Chebyshev series.
#
# Parameters:
# f - Function to fit, expressed as a command prefix to which
# the abscissa will be appended. The function is expected
# to return the ordinate.
# a,b - Interval over which to fit the function
# epsilon - Desired absolute error bound of the fit (approximate)
# n - Maximum order of the Chebyshev polynomial to include.
#
# Results:
# Returns a Tcl command prefix that, when evaluated with an
# abscissa postpended, will give an approximation to the function.
#
# Side effects:
# Evaluates the given function n times, so has whatever side
# effects it has. Also may cache a table of cosines.
#
#----------------------------------------------------------------------
proc ::math::approximate::cheby_fit { f a b {epsilon 1e-8} { n 50 } } {
set halfwidth [expr { ($b - $a) / 2. }]
set midpoint [expr { ( $b + $a ) / 2. }]
set cs [cos_table $n]
set c {}
for { set k 0 } { $k < $n } { incr k } {
lappend c 0.
}
set twon [expr { $n + $n }]
set fourn [expr { $twon + $twon }]
for { set k 1 } { $k <= $twon } { incr k 2 } {
set x [expr { [lindex $cs $k] * $halfwidth + $midpoint }]
set cmd $f; lappend cmd $x; set y [uplevel 1 $cmd]
set m 0
for { set j 0 } { $j < $n } { incr j } {
lset c $j [expr { [lindex $c $j] + $y * [lindex $cs $m] }]
incr m $k
if { $m >= $fourn } {
set m [expr { $m - $fourn }]
}
}
}
set r {}
foreach d $c {
lappend r [expr { $d * 2. / $n }]
}
lset r 0 [expr { [lindex $r 0] / 2. }]
for { set i [expr { $n-1 }] } { $i > 0 } { incr i -1 } {
if { abs([lindex $r $i]) > $epsilon } break
}
return [list [namespace which -command cheby_eval] [lrange $r 0 $i] $a $b]
}
#----------------------------------------------------------------------
#
# cos_table --
#
# Builds a cosine table for Chebyshev fitting.
#
# Parameters:
# n - Number of points.
#
# Results:
# Returns a list of length 4n. The ith element of the list
# is cos( i * pi / ( 2 * n ) )
#
# Side effects:
# Caches the cosine table.
#
#----------------------------------------------------------------------
proc ::math::approximate::cos_table { n } {
variable costable
if { [info exists costable($n)] } {
return $costable($n)
}
set theta [expr { 0.78539816339744828 / $n }]
set s [expr { sin($theta) }]
set alpha [expr { 2. * $s * $s }]
set beta [expr { sin( 2. * $theta ) }]
set c 1.
set s 0.
set table {}
for { set i 0 } { $i + $i <= $n } { incr i } {
lappend ctable $c
lappend stable $s
foreach { c s } \
[list \
[expr { $c - ( $alpha * $c + $beta * $s ) }] \
[expr { $s - ( $alpha * $s - $beta * $c ) }]] break
}
incr i [expr { ( $n % 2 ) - 1 }]
for { incr i -1 } { $i > 0 } { incr i -1 } {
lappend ctable [lindex $stable $i]
lappend stable [lindex $ctable $i]
}
foreach c $stable {
set c [expr {- $c}]
lappend ctable $c
}
foreach c $ctable {
set c [expr { - $c }]
lappend ctable $c
}
if { $n <= 100 } {
set costable($n) $ctable
}
return $ctable
}
if { ! [info exists ::argv0] || ![string equal [info script] $::argv0] } {
return
}
set tcl_precision 17
proc f { x } {
return [expr { pow( sin( 6.28318 * $x - 1.2 ), 5) }]
}
package require Tk
grid [canvas .c -width 600 -height 400 -bg white]
.c create line 20 20 20 380
.c create line 20 380 580 380
proc cx { x } {
expr { 560. * $x + 20. }
}
proc cy { y } {
expr { 200 - 180. * ( $y ) }
}
for { set xx 0 } { $xx <= 10 } { incr xx } {
set x [expr { $xx / 10. }]
.c create text [cx $x] 382 -anchor n -text [format %.1f $x]
}
for { set yy -10 } { $yy <= 10 } { incr yy 2 } {
set y [expr { $yy / 10. }]
.c create text 18 [cy $y] -anchor e -text [format %.1f $y]
}
set l {.c create line}
for { set x 0. } { $x < 1.0025 } { set x [expr { $x + 0.005 }] } {
set cmd f; lappend cmd $x; set y [eval $cmd]
lappend l [cx $x] [cy $y]
}
lappend l -width 5 -fill gray85
eval $l
.c create text 590 10 -anchor ne -tags n -font {Helvetica 12}
after 1000 set done 1; vwait done
set apx0 [::math::approximate::cheby_fit f 0. 1. 1.e-8 120]
for { set n 0 } { $n < [llength [lindex $apx0 1]] } { incr n } {
.c itemconfigure n -text $n
set apx $apx0
lset apx 1 [lrange [lindex $apx 1] 0 $n]
set l {.c create line}
set maxy2 0.
for { set x 0. } { $x < 1.0025 } { set x [expr { $x + 0.005 }] } {
set cmd $apx; lappend cmd $x; set y [eval $cmd]
lappend l [cx $x] [cy $y]
}
lappend l -tags apx
catch { .c delete apx }
eval $l
after 2000 set done 1; vwait done
}
