[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)}]"

----
[[
[Category Mathematics]
| [Category Numerical Analysis]
]]