[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 } ---- [[ Category Mathematics] ]]