Arjen Markus (10 may 2004) Playing with mathematical expressions, I thought I'd make a small procedure to deal with mathematical limits:
sin(x) lim ------ = 1 x->0 x
As I have never seen a textbook handling that kind of things, here is a home-made numerical method. Not tested extensively and the parameter values are not the result of a thorough analysis - quite the opposite: just figured they would be nice values ...
# limit.tcl --
# Tentative method to determine the limit of some mathematical
# expression with one variable nearing a limit value
#
# Note:
# This is very much a do-it-yourself method - I have
# never seen a treatise on the subject of numerically
# calculating limits.
#
# General idea:
# We can use this in the script for evaluating mathematical
# expressions:
# math {
# y = lim x->0 sin(x)/x
# z = lim x->inf (sqrt(x+2)-sqrt(x+1))/(sqrt(x+1)-sqrt(x))
# }
#
# Other "special" forms:
# y = deriv x=0 log(1+x*x)
# y = int x=0 to 1 sin(x)
# y = sum x=1 to 10 1.0/pow(x,2)
#
namespace eval ::math::calculus {
namespace export limit
namespace eval v {} ;# Private namespace
}
# limit --
# Determine the limit of an expression as the variable nears
# the given limit
#
# Arguments:
# varname Name of the variable to reach the limit
# limval The limit value ("inf" is infinite)
# expression The expression to evaluate (make sure it
# is enclosed in braces!)
# level Optional (private use only): number of
# levels to go up
# Result:
# Estimated limit
# Note:
# The expression is evaluated in the caller's scope,
# so that a call like:
# limit x 0 {sin($a*$x)/$x}
# with "a" as a "free" parameter in the calling procedure
#
# Todo:
# It is not yet dealing with the possibility that the
# numerical evaluation may fail due to the large value
# of an argument (say: exp(1.0e30))
#
proc ::math::calculus::limit {varname limval expression {level 1}} {
#
# Prepare the expression
#
set expression [string map [list \$$varname \$::math::calculus::v::$varname] $expression]
#
# It is easiest to just evaluate the expression
# in the caller's scope and do all logic here
#
set dvar 0.25
set factor 0.5
set base 1.0
set limv2 $limval
if { $limval == 0.0 } {
set limv2 1.0
set base 0.0
}
if { $limval == "inf" } {
set limv2 1.0 ;# start carefully
set factor 16.0
}
if { $limval == "-inf" } {
set limv2 -1.0 ;# start carefully
set factor 16.0
}
set count 0
set fvalues {}
set vvalues {}
while { $count < 10 } {
set v::$varname [expr {$limv2*($base+$dvar)}]
if { [catch {
lappend fvalues [uplevel $level [list expr $expression]]
lappend vvalues [set v::$varname]
}] } {
break
}
incr count
set dvar [expr {$dvar*$factor}]
}
#
# Next step:
# Estimate the limiting value via extrapolation
# (use a quadratic polynomial: y = a + b x + c x**2)
#
foreach {x1 x2 x3} [lrange $vvalues end-2 end] {break}
foreach {y1 y2 y3} [lrange $fvalues end-2 end] {break}
if { $limval == "inf" || $limval == "-inf" } {
set x1 [expr {1.0/$x1}]
set x2 [expr {1.0/$x2}]
set x3 [expr {1.0/$x3}]
set limval 0.0
}
set yy [expr {($y1-$y2)*($x1-$x3)-($y1-$y3)*($x1-$x2)}]
set c [expr {$yy /(($x1-$x2)*($x1-$x3)*($x2-$x3))}]
set b [expr {(($y1-$y2)-$c*($x1*$x1-$x2*$x2))/($x1-$x2)}]
set a [expr {$y1-$b*$x1-$c*$x1*$x1}]
return [expr {$a+$b*$limval+$c*$limval*$limval}]
}
#
# Test
#
namespace import ::math::calculus::*
set a 2.0
puts "[limit x 0 {sin($a*$x)/$x}] (expected: $a)"
puts "[limit x inf {atan($x)}] (expected: [expr {acos(-1)/2.0}])"
puts "[limit x 0 {log($x)*$x}] (expected: 0.0)"
puts "[limit x 1 {log($x)/($x-1)}] (expected: 1.0)"
set b 10
puts "[limit x inf {(sqrt($x+$b)-sqrt($x+1))/(sqrt($x+1)-sqrt($x))}] (expected: [expr {$b-1.0}])"
puts "Likely failure:"
puts "[limit x inf {exp($x)/(exp($x)-1)}]"]