a**2 + 2*a*b + b**2 ==> (a+b)**2One kind of manipulation is to find the derivative of some function:
f(x) = x**2 + 2*x df --(x) = 2*x + 2 dxThe rules for that kind of manipulation are simple and can applied mechanically (in contrast: finding the primitive is a challenge that quite often leads to completely new classes of functions!)A very nice technique is to let the expression do the hard work: this is called automatic differentiation and here is one way to do it in Tcl:
- First we define a bunch of procedures that allow us to express the function f above as a nesting of ordinary Tcl commands (these are easier to manipulate than the original expression):
proc f {x} {
return [+ [* $x $x] [* 2 $x]]
}- Then we use a list of two elements instead of plain numbers as arguments:
set x [list 2.0 1.0]
- The meaning of the first is: the value of x, the meaning of the second is: the value of the derivative of x with respect to the independent variable x.
- Our procedures manipulate these two numbers, instead of the single value - this leads to the automatic evaluation of the function and its derivative. (Constants are represented as single numbers or as lists with 0.0 as the second element).
f(x) = g(x) + h(x) ==> f'(x) = g'(x) + h'(x) f(x) = g(x) * h(x) ==> f'(x) = g'(x)*h(x) + g(x)*h'(x)etc. can be computed without us having to do more than specifying the way operators like + and * behave under differentiation.
namespace eval ::derivs {
proc + {x y} {
if { [llength $x] == 1 } {
set x [list $x 0.0]
}
if { [llength $y] == 1 } {
set y [list $y 0.0]
}
set v [expr {[lindex $x 0]+[lindex $y 0]}]
set dv [expr {[lindex $x 1]+[lindex $y 1]}]
return [list $v $dv]
}
proc * {x y} {
if { [llength $x] == 1 } {
set x [list $x 0.0]
}
if { [llength $y] == 1 } {
set y [list $y 0.0]
}
set v [expr {[lindex $x 0]*[lindex $y 0]}]
set dv [expr {[lindex $x 1]*[lindex $y 0] +
[lindex $x 0]*[lindex $y 1]}]
return [list $v $dv]
}
# - and / are left as an exercise ...
#
# Also define a few functions to make it more interesting
#
foreach {f df} {sin cos exp exp cos -sin log 1.0/ sqrt 0.5/sqrt} {
proc $f {x} [string map [list F $f DF $df] {
if { [llength $x] == 1 } {
set x [list $x 0.0]
}
set v [expr {F([lindex $x 0])}]
set dv [expr {[lindex $x 1]*DF([lindex $x 0])}]
return [list $v $dv]
}]
}
}
#
# A few tests:
# Define two (global) functions
#
proc f {x} {
return [* $x [+ 2 $x]]
}
proc g {x} {
return [* 2 [* [sin $x] [cos $x]]]
}
namespace eval ::derivs {
#
# f'(x) = (2+x)+x = 2+2x
# g'(x) = 2*cos(x)**2-2*sin(x)**2 = 2*cos(2x)
#
# "Import" the two functions - well, this is one way to do it :)
#
foreach p {f g} {
proc $p [info args ::$p] [info body ::$p]
}
#
# Test the derivatives
#
puts "Parabola:"
foreach x {0.0 0.5 1.0 1.5 2.0} {
set xl [list $x 1.0] ;# Very important!!
puts "[f $xl] -- [expr {2.0+2.0*$x}]"
}
puts "sin(2x):"
foreach x {0.0 0.5 1.0 1.5 2.0} {
set xl [list $x 1.0]
puts "[g $xl] -- [expr {2.0*cos(2.0*$x)}]"
}
}This can be extended to the evaluation of second derivatives if you like, or you can construct the expression that evaluates the derivative.Another, somewhat simpler, variation is to determine the size of the expression tree - something that is important with genetic programming techniques. Merely evaluate the above functions in a different namespace ... a neat trick Donal Fellows mentioned on comp.lang.tcl: namespace eval ::treesize {
foreach op {+ * - /} {
#
# Only the "arity" counts ...
#
proc $op {x y} {
set count 1
if { [llength $x] == 1 } {
incr count
} else {
incr count [lindex $x 1]
}
if { [llength $y] == 1 } {
incr count
} else {
incr count [lindex $y 1]
}
return [list COUNT $count]
}
}
#
# Also define a few functions to make it more interesting
#
foreach f {sin exp cos} {
proc $f {x} {
set count 1
if { [llength $x] == 1 } {
incr count
} else {
incr count [lindex $x 1]
}
return [list COUNT $count]
}
}
}
#
# Reuse the functions above ...
#
namespace eval ::treesize {
#
# "Import" the two functions
#
foreach p {f g} {
proc $p [info args ::$p] [info body ::$p]
}
puts "Size of 'f': [f 1]"
puts "Size of 'g': [g 1]"
}So, with the proper operators and mathematical functions we need only to define the actual expression once to automatically derive all kinds of mathematical and symbolic properties!Suppose you have an expression like:f(x) = ln( x + sqrt(1+x**2) )and you need to compute the derivative, that is trivial with the code above. Manually verifying that its first derivative is:
1
f'(x) = ------------
sqrt(1+x**2)and then computing the value would take some time :)TV From an earlier look at this page I recall the algorithm has limitation with respect to the types of forms of equations? If it is a general derivative routine, it would be good, regardless of the formula format.Don't forget Maxima can also do derivations symbolically, quite powerfully, and even the reverse, which is a lot harder still, and it can be driven by tcl.AM The method presented above is quite general, the implementation may be a bit limited due to my laziness :). The very attractive thing about it is that you do not need to resort to numerical differentiation which suffers from accuracy problems, not to mention problems with singularities ... This is an exact method with no caveats beyond the mathematical - I mean: the derivative has to exist.
SYMDIFF
is a tool for symbolic differentiation. Features:- Expression parser
- Symbolic differentiation of expressions with respect to arbitrary number of variables.
- User defined differentiation rules for arbitrary functions
- Common subexpression elimination for group of expressions
- C++ library
- Tcl scripting interface
