Symbolic differentiation of a formula, then numeric evaluation.

The program roots back to I believe work of D.H.D. Warren and colleagues in the late 1970s; variants are part of the Prolog culture.

This program shows the transformation of an engineered language (e.g. trigonometric functions) according to math textbook rules, then handing over the transformed formula to further processing, here numeric evaluation.

It serves as example for the combination of symbolic and numeric solvers, and is the classic example for the power of a recursive program (the differentiation rules) working on a recursive data structure (the math expression).

Here is the source for the differentiation rules:

% diff.min	- symbolic differentiation rules
% diff(+Formula,+AfterVariable,-DifferentiatedFormula)

diff(X1,X2,1) :-  var(X1), !, X1 == X2.
diff(C,X,0) :- constant(C,X),!.
diff(-U,X,-DU) :- diff(U,X,DU).
diff(U+V,X,DU+DV) :- diff(U,X,DU), diff(V,X,DV).
diff(U-V,X,DU-DV) :- diff(U,X,DU), diff(V,X,DV).
diff(C*U,X,C*DU) :- constant(C,X), diff(U,X,DU).
diff(U*V,X,DV*U+DU*V) :- \+ constant(U,X),!, diff(U,X,DU), diff(V,X,DV).
diff(U/V,X,DU) :- diff(U*V^ (-1),X,DU).
diff(U^C,X,C*DU*U ^ (C-1)) :- constant(C,X),!,diff(U,X,DU).
diff(ln(U),X,DU*U ^ (-1)) :- diff(U,X,DU).
diff(sin(X1),X,cos(X)) :- X1 == X.
diff(cos(X1),X,-sin(X)) :- X1 == X.
diff(exp(X1),X,exp(X)) :- X1 == X.
diff(GoF,X,DG*DF) :- 
	GoF =.. [G,F],
	known_function(G),
	Y =..[G,U],
	diff(Y,U,DG), 
	diff(F,X,DF),
	U = F.
diff(F,X,_DF) :- throw('Sorry, no rule for diff'(F,X)).

known_function(G) :- member(G,[ln,sin,cos,exp]).

% constant(+C,X) <- C is a constant expression with respect to X
constant(C,_X) :- atomic(C),!.
constant(C,_X) :- compound(C),
	C = const(_NameOfSymbolicConstant),!.	% a declared constant
constant(C,X) :-
	compound(C), C =.. [BinOp,A,B],!,
	binary_operator(BinOp), constant(A,X), constant(B,X).
constant(C,X) :-
	compound(C), C =.. [UnOp,A],!,
	unary_operator(UnOp), constant(A,X).

% binary_operator(+BinOp) <- BinOp is a binary operator
binary_operator(BinOp) :- member(BinOp,[-,+,*,/,^]).

% unary_operator(+UnOp) <- UnOp is a unary operator
unary_operator(UnOp) :-	member(UnOp,[-,ln,exp,sin,cos]).