Example: A Symbolic Differentiator - Maple Programming Help

Home : Support : Online Help : Applications and Example Worksheets : Language and System : examples/SymbolicDifferentiator

Example: A Symbolic Differentiator

Introduction

 This section illustrates the various module concepts through a symbolic differentiator example. Since Maple provides a built-in differentiator diff, the example symbolic differentiator is named differentiate. Its (final) implementation is in the module DiffImpl (later in this chapter), which holds all the local state for the program. Much of the code for the differentiator is designed to implement either a standard rule (such as the rule that the derivative of a sum is the sum of the derivatives of the summands), or special case rules for mathematical functions such as sin and exp. The example differentiator handles only real valued functions of a single real variable.
 The following example shows several steps in the development of the module, from a very simple first try to the final, fully functional program. The final form of the differentiator is a good illustration of a very common Maple design pattern. This pattern arises when you have a single top-level routine that dispatches a number of subroutines to handle special cases using special purpose algorithms.

The First Attempt

 This initial example presents the differentiator as an ordinary procedure, not a module.
 > differentiate := proc( expr, var )     local a, b;     if type( expr, 'constant' ) then         0;     elif expr = var then         1;     elif type( expr, '+' ) then         map( procname, args );     elif type( expr, '^' ) then         a, b := op( expr );         if a = var and not has( b, var ) then             b * a ^ ( b - 1 );         else             'procname( args )';         end if;     elif type( expr, '*' ) then         a, b := op( 1, expr ), subsop( 1 = 1, expr );         procname( a, var ) * b + a * procname( b, var );     else         'procname( args )';     end if; end proc:
 Trivial cases are handled first: The derivative of a constant expression is equal to 0, and the derivative of the variable with respect to which we are differentiating is equal to 1. The additivity of the derivative operator is expressed by mapping the procedure over sums, using the command
 > map( procname, args );
 ${\mathrm{procname}}{}\left({\mathrm{args}}\right)$ (1)
 This is commonly used to map a procedure over its first argument, passing along all the remaining arguments. Only the simple case of powers of the differentiation variable is handled so far, provided also that the power is independent of the differentiation variable. The product rule for derivatives is expressed by splitting expressions of type product into two pieces:
 • the first factor in the product, and
 • the product of all the remaining factors.
 This is achieved by the double assignment of
 > a, b := op( 1, expr ), subsop( 1 = 1, expr );
 ${a}{,}{b}{≔}{\mathrm{expr}}{,}{1}$ (2)
 so the input expression expr is expressed as expr = a * b. The standard technique of returning unevaluated is used so that computation can proceed symbolically on expressions that the procedure is unable to differentiate.
 This first example is simple, but it is already able to handle polynomials with numeric coefficients.
 > differentiate( 2 - x + x^2 + 3*x^9, x );
 ${27}{}{{x}}^{{8}}{+}{2}{}{x}{-}{1}$ (3)
 However, it fails on expressions containing calls to standard mathematical functions.
 > differentiate( sin( x ), x );
 ${\mathrm{differentiate}}{}\left({\mathrm{sin}}{}\left({x}\right){,}{x}\right)$ (4)
 It is also unable to deal successfully with symbolic coefficients.
 > differentiate( a*x^2 + b*x + c, x );
 ${\mathrm{differentiate}}{}\left({\mathrm{expr}}{,}{x}\right){}{{x}}^{{2}}{+}{2}{}{\mathrm{expr}}{}{x}{+}{\mathrm{differentiate}}{}\left({c}{,}{x}\right){+}{1}$ (5)

 To add the missing functionality, add a case for expressions of type function.
 > differentiate := proc( expr, var )     local a, b;     if not has( expr, var ) then         0;     elif expr = var then         1;     elif type( expr, '+' ) then         map( procname, args );     elif type( expr, '^' ) then         a, b := op( expr );         if not has( b, var ) then             b * a ^ ( b - 1 ) * procname( a, var );         else             'procname( args )';         end if     elif type( expr, '*' ) then         a, b := op( 1, expr ), subsop( 1 = 1, expr );         procname( a, var ) * b + a * procname( b, var )     elif type( expr, 'function' ) and nops( expr ) = 1 then         # functions of a single variable; chain rule         b := op( 0, expr ); # the name of the function         a := op( 1, expr ); # the argument         if b = 'sin' then             cos( a ) * procname( a, var );         elif b = 'cos' then             -sin( a ) * procname( a, var );         elif b = 'exp' then             exp( a ) * procname( a, var );         elif b = 'ln' then             ( 1 / a ) * procname( a, var );         else             'procname( args )';         end if;     else         'procname( args )';     end if; end proc:
 This uses the chain rule to compute the derivatives of calls to known functions.
 > differentiate( sin( x ) + cos( exp( x ) ), x );
 ${\mathrm{cos}}{}\left({x}\right){-}{\mathrm{sin}}{}\left({{ⅇ}}^{{x}}\right){}{{ⅇ}}^{{x}}$ (6)
 > differentiate( sin( x^2 ) + cos( x^2 ), x );
 ${2}{}{\mathrm{cos}}{}\left({{x}}^{{2}}\right){}{x}{-}{2}{}{\mathrm{sin}}{}\left({{x}}^{{2}}\right){}{x}$ (7)
 > differentiate( sin( x )^2 + cos( x )^3, x );
 ${2}{}{\mathrm{sin}}{}\left({x}\right){}{\mathrm{cos}}{}\left({x}\right){-}{3}{}{{\mathrm{cos}}{}\left({x}\right)}^{{2}}{}{\mathrm{sin}}{}\left({x}\right)$ (8)
 At the same time, this has also improved the handling of expressions independent of the variable of differentiation.
 > differentiate( a*x^2 + b*x + c, x );
 ${2}{}{\mathrm{expr}}{}{x}{+}{1}$ (9)
 This is effected by using the expression has( expr, var ) instead of the weaker test type( expr, 'constant' ). The power rule now handles more than just powers of var.
 > differentiate( sin( x )^2, x );
 ${2}{}{\mathrm{sin}}{}\left({x}\right){}{\mathrm{cos}}{}\left({x}\right)$ (10)
 However, adding new functions to the differentiator is tedious and error prone, and the job of handling the chain rule must be repeated for each function recognized by it.

Introducing a Function Table

 Many functions (that you need to add) and the rules used for their differentiation can be stored in a table as follows:
 > differentiate := proc( expr, var )     local a, b, functab;     functab := table();     functab[ 'sin' ] := 'cos';     functab[ 'cos' ] := x -> -sin( x );     functab[ 'exp' ] := exp;     functab[ 'ln' ] := x -> 1 / x;     if not has( expr, var ) then         0;     elif expr = var then         1;     elif type( expr, '+' ) then         map( procname, args );     elif type( expr, '^' ) then         a, b := op( expr );         if a = var and not has( b, var ) then             b * a ^ ( b - 1 ) * procname( a, var );         else             'procname( args )';         end if     elif type( expr, '*' ) then         a, b := op( 1, expr ), subsop( 1 = 1, expr );         procname( a, var ) * b + a * procname( b, var );     elif type( expr, 'function' ) and nops( expr ) = 1 then         # functions of a single variable; chain rule         b := op( 0, expr ); # the name of the function         a := op( 1, expr ); # the argument         if assigned( functab[ b ] ) then             # This is a \"known\" function             functab[ b ]( a ) * procname( a, var );         else             # This function is not known; return unevaluated             'procname( args )';         end if;     else         'procname( args )';     end if; end proc:
 This not only simplifies the code used for the function case, but also makes it very easy to add new functions.

Drawbacks

 Unfortunately, this implementation has serious drawbacks.
 • It is not extensible. The known functions are hard-coded as part of the procedure definition for differentiate. New functions cannot be added without editing this source code.
 • A second problem relates to performance. A complete implementation would require a table of dozens or hundreds of functions. That large table would need to be created and initialized each time differentiate is invoked.

Encapsulation and Extensibility

 One way to fix both problems is to make the table of functions a global variable. However, using global variables can be dangerous, because they pollute the user namespace and are subject to unwanted inspection and tampering.

Solution

 A better solution is to put the differentiate procedure, along with its table of functions, into a module. The table is then initialized only once--when the module is created--and can be saved to a Maple repository with the rest of the module by using a savelib call. By making the table a local variable of the module, you prevent users from modifying the table or otherwise inspecting it in unwanted ways.
 This does not prevent you from making the differentiator user-extensible, however. You can add an access procedure addFunc that allows users to add rules for differentiating new functions. For example, you can use the call
 ${\mathrm{addFunc}}{}\left({\mathrm{cos}}{,}{x}{↦}{-}{\mathrm{sin}}{}\left({x}\right)\right)$ (11)
 to add the derivative of the sin function. The export addFunc of the DiffImpl module is a procedure that requires two arguments. The first is the name of a function whose derivative is to be made known to the differentiator. The second is a Maple procedure of one argument that expresses the derivative of the function being added.
 With this strategy in mind, you can create a module DiffImpl, with principal export differentiate. At the same time, you can also make the basic differentiation rules extensible.
 Here is the complete source code for the differentiator with these improvements.
 > DiffImpl := module()     description "a symbolic differentiator";     local    functab, ruletab, diffPower;     export    differentiate, addFunc, addRule, rule;     addFunc := proc( fname::symbol, impl )         functab[ fname ] := impl;     end proc;     addRule := proc( T, impl )         if type( T, '{ set, list }' ) then             map( procname, args );         elif type( T, 'And( name, type )' ) then             ruletab[ T ] := impl;         else             error "expecting a type name, but got %1", T;         end if     end proc;     rule := proc( T )         if type( T, 'And( name, type )' ) then             if assigned( ruletab[ T ] ) then                 eval( ruletab[ T ], 1 );             else                 error "no rule for expressions of type %1", T;             end if;         else             error "expecting a type symbol, but got %1", T;         end if;     end proc;     differentiate := proc( expr, var )         local a, b, e;         if not has( expr, var ) then             0;         elif expr = var then             1;         elif type( expr, 'function' ) and nops( expr ) = 1 then             e := op( 0, expr );             a := op( expr );             if assigned( functab[ e ] ) then                 functab[ e ]( a ) * procname( a, var );             else                 'procname( args )';             end if;         else             b := whattype( expr );             if assigned( ruletab[ b ] ) then                 ruletab[ b ]( expr, var );             else                 'procname( args )';             end if;         end if;     end proc;     addRule( '{list,set,tabular}',             () -> map( differentiate, args ) );     addRule( '+',              () -> map( differentiate, args ) );     addRule( '*',       (expr,var) ->         op(1,expr)*differentiate(subsop(1=1,expr),var)            + differentiate(op(1,expr),var)*subsop(1=1,expr) );     diffPower := proc( expr, var )         local    b, e;         Assert( type( expr, '^' ) );         b, e := op( expr );         if has( e, var ) then             expr * ( differentiate( e, var ) * ln( b )                + e * differentiate( b, var ) / b );         else # simpler formula             e * b^(e - 1) * differentiate( b, var );         end if;     end proc;     addRule( '^', eval( diffPower ) );     addFunc( 'sin', cos );     addFunc( 'cos', x -> -sin(x) );     addFunc( 'exp', exp );     addFunc( 'ln', x -> 1/x );     # ... etc. end module:
 > differentiate := DiffImpl:-differentiate:
 To give the set of rules for nonfunctional expressions similar extensibility, you can store those rules in a table. The table is indexed by the primary (or basic) type name for the expression type, as given by the Maple procedure whattype.
 > whattype( a + 2 );
 ${\mathrm{+}}$ (12)
 > whattype( a / b );
 ${\mathrm{symbol}}$ (13)
 > whattype( a^sqrt(2) );
 ${\mathrm{^}}$ (14)
 > whattype( [ f( x ), g( x ) ] );
 ${\mathrm{list}}$ (15)
 A rule is expressed by a procedure of two arguments, expr and var, in which expr is the expression to be differentiated, and var is the variable of differentiation. For instance, to make the differentiator handle items such as sets and lists by differentiating their individual components, add the rule
 > addRule('{list, set, tabular}', () -> map( differentiate, args) );
 ${\mathrm{addRule}}{}\left(\left\{{\mathrm{list}}{,}{\mathrm{set}}{,}{\mathrm{tabular}}\right\}{,}\left(\right){↦}{\mathrm{map}}{}\left({\mathrm{differentiate}}{,}{\mathrm{args}}\right)\right)$ (16)
 The first version of the differentiator dealt with sums by mapping itself over the sum expression. In the new scheme, this is expressed by the statement in the module body. The advantage of using this scheme is that users, as well as the author of the differentiator, can extend the system. Having instantiated the module DiffImpl, any user can add rules or new functions, simply by issuing appropriate calls to addRule and addFunc.
 The differentiator cannot handle the procedure tan.
 > differentiate( tan( x )/exp( x ), x );
 ${-}\frac{{\mathrm{tan}}{}\left({x}\right)}{{{ⅇ}}^{{x}}}{+}\frac{{\mathrm{differentiate}}{}\left({\mathrm{tan}}{}\left({x}\right){,}{x}\right)}{{{ⅇ}}^{{x}}}$ (17)
 You must add it to the database of known functions.
 > DiffImpl:-addFunc( 'tan', x -> 1 + tan(x)^2 );
 ${x}{↦}{1}{+}{{\mathrm{tan}}{}\left({x}\right)}^{{2}}$ (18)
 > differentiate( tan( x )/exp( x ), x );
 ${-}\frac{{\mathrm{tan}}{}\left({x}\right)}{{{ⅇ}}^{{x}}}{+}\frac{{1}{+}{{\mathrm{tan}}{}\left({x}\right)}^{{2}}}{{{ⅇ}}^{{x}}}$ (19)
 Similarly, there is not yet any rule for handling equations and other relations.
 > differentiate( y( x ) = sin( x^2 ) - cos( x^3 ), x );
 ${\mathrm{differentiate}}{}\left({y}{}\left({x}\right){=}{\mathrm{sin}}{}\left({{x}}^{{2}}\right){-}{\mathrm{cos}}{}\left({{x}}^{{3}}\right){,}{x}\right)$ (20)
 > DiffImpl:-addRule( '{ =, <, <= }',               () -> map( differentiate, args ) );
 $\left\{\left(\right){↦}{\mathrm{map}}{}\left({\mathrm{differentiate}}{,}{\mathrm{args}}\right)\right\}$ (21)
 > differentiate( y( x ) = sin( x^2 ) - cos( x^3 ), x );
 ${\mathrm{differentiate}}{}\left({y}{}\left({x}\right){,}{x}\right){=}{2}{}{\mathrm{cos}}{}\left({{x}}^{{2}}\right){}{x}{+}{3}{}{\mathrm{sin}}{}\left({{x}}^{{3}}\right){}{{x}}^{{2}}$ (22)

The Extension Mechanism is Module Aware

 Do not confuse the extension mechanism previously proposed for the differentiator with the extension mechanism used by the built-in Maple command diff . The diff command uses a traditional string concatenation mechanism for adding knowledge of the derivatives of functions, and all its rules are built-in, so they cannot be extended. For instance, to add a new function F to the Maple built-in diff command, you can define a procedure diff/F that computes the derivative of F.
 By contrast, the extension mechanism used in the differentiate example is module aware. To add knowledge of the derivative of some top-level function F, you can issue a command, such as
 > DiffImpl:-addFunc( 'F', x -> sin( x ) + cos( x ) );
 ${x}{↦}{\mathrm{sin}}{}\left({x}\right){+}{\mathrm{cos}}{}\left({x}\right)$ (23)
 The derivative of F( x ) is sin( x ) + cos( x ).) Define a module with some special functions, one of which is also called F.
 > SpecFuncs := module()     export F; # etc.     # definition of F() and others end module:
 You can now add this new F to the known functions.
 > DiffImpl:-addFunc( SpecFuncs:-F, x -> exp( 2 * x ) );
 ${x}{↦}{{ⅇ}}^{{2}{\cdot }{x}}$ (24)
 > differentiate( F( x ), x );
 ${\mathrm{sin}}{}\left({x}\right){+}{\mathrm{cos}}{}\left({x}\right)$ (25)
 > use SpecFuncs in     differentiate( F( x ), x ); end use;
 ${{ⅇ}}^{{2}{}{x}}$ (26)
 With the traditional mechanism, this does not work.
 > cat( diff/, F ) := x -> sin( x ) + cos( x );
 ${\mathrm{cat}}{}\left({\mathrm{diff/}}{,}{F}\right){≔}{x}{↦}{\mathrm{sin}}{}\left({x}\right){+}{\mathrm{cos}}{}\left({x}\right)$ (27)
 > diff( F( x ), x );
 ${\mathrm{sin}}{}\left({x}\right){+}{\mathrm{cos}}{}\left({x}\right)$ (28)
 > use SpecFuncs in     cat( diff/, F ) := x -> exp( 2 * x );     diff( F( x ), x ); end use;
 ${{ⅇ}}^{{2}{}{x}}$ (29)
 The definition for the global F has been lost.
 > diff( F( 2 * x ), x );
 ${{ⅇ}}^{{4}{}{x}}$ (30)
 (You can use a different argument to diff to avoid recalling the answer from its remember table.) The traditional mechanism fails because it relies on the external representation of names, and not upon their bindings, so each attempt to define an extension to diff in fact adds a definition for the derivative of all functions whose names are spelled "F".
 Note: A commented version of the differentiate module is available in the samples/ProgrammingGuide directory of the Maple installation. The implementation shown in the text has been somewhat simplified.