Symbolic Differentiation 

The first step in symbolic algebra is defining a representation. Getting expressions into a suitable form is actually straightforward; no parsing of expressions is needed, since we have Lua to do that for us. Using the pl.func
library does all the hard work; it redefines the arithmetic operations to work on placeholder expressions (PEs), which are Lua expressions involving dummy variables called placeholders. pl.func
defines standard placeholders for arguments called _1
,_2
, etc but the Var
function will create new ones of our chosing:
utils.import 'pl.func' a,b,c,d = Var 'a,b,c,d' print(a+b+c+d)
Which will indeed print out the expression in a readable form. PE operator expressions are stored as combinations of tables looking like {op='+',x,y}
, which have an associated metatable which defines the metamethods like __add
and so forth. As a tree, with the usual associativity of Lua operators we get:
It is irritating to draw these diagrams, so a better notation is Lispstyle Sexpressions:
1: (+ (+ (+ a b) c) d)
However, with the various manipulations we will perform, this canonical form is not the only possible representation of a+b+c+d
:
2: (+ a (+ b (+ c d))) 3: (+ (+ a b) (+ c d))
Now, experience shows that this leads to madness. Instead, it's easier to go for the canonical Lisp representation:
4: (+ a b c d)
Many operations become straightforward once this is in place, for instance comparing with (+ a c b d)
is just a matter of doing a 'compare with no order' on the arguments. Displaying PEs in this form is straightforward. isPE
simply checks the expression to see if it is a placeholder expression, by looking at the metatable. PEs with op=='X'
are placeholder variables, so the rest must be expression nodes.
function sexpr (e) if isPE(e) then if e.op ~= 'X' then local args = tablex.imap(sexpr,e) return '('..e.op..' '..table.concat(args,' ')..')' else return e.repr end else return tostring(e) end end
The first task is to balance the expressions, which converts representations 13 into 4.
function balance (e) if isPE(e) and e.op ~= 'X' then local op,args = e.op if op == '+' or op == '*' then args = rcollect(e) else args = imap(balance,e) end for i = 1,#args do e[i] = args[i] end end return e end
For the noncommutative operators, the idea is just to balance all the subexpressions by mapping balance
over the array part of the PE, which is the argument list. These are then copied back inplace. The nontrivial part is dealing with + and *, where it is necessary to collect all the arguments from expression trees looking like 1,2 or 3 and convert them into the fourth form.
function tcollect (op,e,ls) if isPE(e) and e.op == op then for i = 1,#e do tcollect(op,e[i],ls) end else ls:append(e) return end end function rcollect (e) local res = List() tcollect(e.op,e,res) return res end
This recursively goes down sameoperator chains (the (+ (+ ...)
mentioned earlier) and collects the arguments, flattening them into nary + or * expressions.
Here is a useful function, which follows the same recursive pattern:
 does this PE contain a reference to x? function references (e,x) if isPE(e) then if e.op == 'X' then return x.repr == e.repr else return find_if(e,references,x) end else return false end end
Here are functions to create nary products and sums:
function muli (args) return PE{op='*',unpack(args)} end function addi (args) return PE{op='+',unpack(args)} end
With this in place, the basic differentiation rules are not difficult. Firstly, only consider subexpressions which do contain the variable:
function diff (e,x) if isPE(e) and references(e,x) then local op = e.op if op == 'X' then return 1 else local a,b = e[1],e[2] if op == '+' then  differentiation is linear local args = imap(diff,e,x) return balance(addi(args)) elseif op == '*' then  product rule local res,d,ee = {} for i = 1,#e do d = fold(diff(e[i],x)) if d ~= 0 then ee = {unpack(e)}  make a copy ee[i] = d append(res,balance(muli(ee))) end end if #res > 1 then return addi(res) else return res[1] end elseif op == '^' and isnumber(b) then  power rule return b*x^(b1) end end else return 0 end end
The derivative of a sum of expressions is the sum of the derivatives. Again, imap
does the job of applying the function recursively over the subexpressions. After constructing the result, we rebalance for luck.
The product rule is given here in its general form, with an explicit check for terms which work out to zero  that is the job of fold
, which will be discussed next.
(uvw..)' = u'vw.. + uv'w... + uvw'... + ...
And finally, the simple power rule. Note how the result can be expressed in a straightforward fashion, since all these operators are acting on PEs.
In fact, all these rules are certainly clearer if you use form 1, binary + and *! But then simplification becomes unbearable. And simplification ('folding') is the tricky one to get right. fold
is a longish function, so I will deal with it in sections:
local op = e.op local addmul = op == '*' or op == '+'  first fold all arguments local args = imap(fold,e) if not addmul and not find_if(args,isPE) then  no placeholders in these args, we can fold the expression. local opfn = optable[op] if opfn then return opfn(unpack(args)) else return '?' end elseif addmul then
The first if
is looking for a case where a subexpression has no symbols, i.e. it is something like 2*5
or 10^2
; in this case, the constant can be completely folded. optable
(defined in pl.operator
) gives a mapping between the operator names and the actual function implementing them.
elseif op == '^' then if args[2] == 1 then return args[1] end  identity if args[2] == 0 then return 1 end end return PE{op=op,unpack(args)}
This clause is clearing up expressions like x^1
and y^0
which naturally arise from the power rule in diff
. Once args
has been processed, the expression can be put together again.
The bulk of this routine handles the awkward twins, + and *.
 split the args into two classes, PE args and nonPE args. local classes = List.partition(args,isPE) local pe,npe = classes[true],classes[false]
List.partition
takes a list and a function, which takes one argument and returns a single value. The result is table where the keys are the returned values, and the values are lists of those elements where the function returned that value. So:
List{1,2,3,4}:partition(function(x) return x > 2 end) > {false={1,2},true={3,4}} List{'one',math.sin,10,20,{1,2}}:partition(type) > {function={function: 00369110},string={one},number={10,20},table={{{1,2}} }
(Mathematically, these are referred to as equivalence classes and partition
would be called the quotient set)
In this case, we want to separate the nonsymbolic arguments from the symbolic arguments; order does not matter. The nonsymbolic arguments npe
can folded into a constant. At this point, operator identity rules can kick in, so that we can drop (* 0 x)
and simplify (+ 0 x)
to be just x
.
The final simplification is replacing repeated values, so that (* x x)
should become (^ x 2)
and (+ x x x)
should become (* x 3)
. count_map
from pl.tablex
will do the job. It is given a listlike table and a function which defines equivalence, and returns a map from the values to the number of their occurrences, so that count_map{'a','b','a'}
is {a=2,b=1}
.
Given this test function:
function testdiff (e) balance(e) e = diff(e,x) balance(e) print('+ ',e) e = fold(e) print(' ',e) end
and these cases:
testdiff(x^2+1) testdiff(3*x^2) testdiff(x^2 + 2*x^3) testdiff(x^2 + 2*a*x^3 + x^4) testdiff(2*a*x^3) testdiff(x*x*x)
we get these results, showing why something like fold
is so necessary to process the result of diff
.
+ 2 * x ^ 1 + 0  2 * x + 3 * 2 * x  6 * x + 2 * x ^ 1 + 2 * 3 * x ^ 2  2 * x + 6 * x ^ 2 + 2 * x ^ 1 + 2 * a * 3 * x ^ 2 + 4 * x ^ 3  6 * a * x ^ 2 + 4 * x ^ 3 + 2 * x + 2 * a * 3 * x ^ 2  6 * a * x ^ 2 + 1 * x * x + x * 1 * x + x * x * 1  x ^ 2 * 3
https://github.com/stevedonovan/Penlight/blob/master/examples/symbols.lua
https://github.com/stevedonovan/Penlight/blob/master/examples/testsymbols.lua