Example: A Symbolic Differentiator - Maple Programming Help

Online Help

All Products    Maple    MapleSim


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 );

procnameargs

(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,bexpr,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 );

27x8+2x1

(3)
  

However, it fails on expressions containing calls to standard mathematical functions.

differentiate( sin( x ), x );

differentiatesinx,x

(4)
  

It is also unable to deal successfully with symbolic coefficients.

differentiate( a*x^2 + b*x + c, x );

differentiateexpr,xx2+2exprx+differentiatec,x+1

(5)

Adding Missing Functionality

  

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 );

cosxsinⅇxⅇx

(6)

differentiate( sin( x^2 ) + cos( x^2 ), x );

2cosx2x2sinx2x

(7)

differentiate( sin( x )^2 + cos( x )^3, x );

2sinxcosx3cosx2sinx

(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 );

2exprx+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 );

2sinxcosx

(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

addFunc('cos', x-> -sin(x));

addFunccos,xsinx

(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 );

`+`

(12)

whattype( a / b );

symbol

(13)

whattype( a^sqrt(2) );

`^`

(14)

whattype( [ f( x ), g( x ) ] );

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) );

addRulelist,set,tabular,mapdifferentiate,args

(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 );

tanxⅇx+differentiatetanx,xⅇx

(17)
  

You must add it to the database of known functions.

DiffImpl:-addFunc( 'tan', x -> 1 + tan(x)^2 );

x1+tanx2

(18)

differentiate( tan( x )/exp( x ), x );

tanxⅇx+1+tanx2ⅇ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 );

differentiateyx=sinx2cosx3,x

(20)

DiffImpl:-addRule( '{ `=`, `<`, `<=` }',
              () -> map( differentiate, args ) );

mapdifferentiate&comma;args

(21)

differentiate( y( x ) = sin( x^2 ) - cos( x^3 ), x );

differentiateyx&comma;x=2cosx2x+3sinx3x2

(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 ) );

xsinx+cosx

(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&ExponentialE;2x

(24)

differentiate( F( x ), x );

sinx+cosx

(25)

use SpecFuncs in
    differentiate( F( x ), x );
end use;

&ExponentialE;2x

(26)
  

With the traditional mechanism, this does not work.

cat( `diff/`, F ) := x -> sin( x ) + cos( x );

catdiff/&comma;Fxsinx+cosx

(27)

diff( F( x ), x );

sinx+cosx

(28)

use SpecFuncs in
    cat( `diff/`, F ) := x -> exp( 2 * x );
    diff( F( x ), x );
end use;

&ExponentialE;2x

(29)
  

The definition for the global F has been lost.

diff( F( 2 * x ), x );

&ExponentialE;4x

(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.

Return to the Index of Example Worksheets.

See Also

module