Math::Symbolic::Derivative - Derive Math::Symbolic trees

Math-Symbolic documentation Contained in the Math-Symbolic distribution.

Index

Code Index:

NAME

Top

Math::Symbolic::Derivative - Derive Math::Symbolic trees

SYNOPSIS

Top

  use Math::Symbolic::Derivative qw/:all/;
  $derived = partial_derivative($term, $variable);
  # or:
  $derived = total_derivative($term, $variable);

DESCRIPTION

Top

This module implements derivatives for Math::Symbolic trees. Derivatives are Math::Symbolic::Operators, but their implementation is drawn from this module because it is significantly more complex than the implementation of most operators.

Derivatives come in two flavours. There are partial- and total derivatives.

Explaining the precise difference between partial- and total derivatives is beyond the scope of this document, but in the context of Math::Symbolic, the difference is simply that partial derivatives just derive in terms of explicit dependency on the differential variable while total derivatives recongnize implicit dependencies from variable signatures.

Partial derivatives are faster, have been tested more thoroughly, and are probably what you want for simpler applications anyway.

EXPORT

None by default. But you may choose to import the total_derivative() and partial_derivative() functions.

CLASS DATA

Top

The package variable %Partial_Rules contains partial derivative rules as key-value pairs of names and subroutines.

SUBROUTINES

Top

partial_derivative

Takes a Math::Symbolic tree and a Math::Symbolic::Variable as argument. third argument is an optional boolean indicating whether or not the tree has to be cloned before being derived. If it is true, the subroutine happily stomps on any code that might rely on any components of the Math::Symbolic tree that was passed to the sub as first argument.

total_derivative

Takes a Math::Symbolic tree and a Math::Symbolic::Variable as argument. third argument is an optional boolean indicating whether or not the tree has to be cloned before being derived. If it is true, the subroutine happily stomps on any code that might rely on any components of the Math::Symbolic tree that was passed to the sub as first argument.

AUTHOR

Top

Please send feedback, bug reports, and support requests to the Math::Symbolic support mailing list: math-symbolic-support at lists dot sourceforge dot net. Please consider letting us know how you use Math::Symbolic. Thank you.

If you're interested in helping with the development or extending the module's functionality, please contact the developers' mailing list: math-symbolic-develop at lists dot sourceforge dot net.

List of contributors:

  Steffen Müller, symbolic-module at steffen-mueller dot net
  Stray Toaster, mwk at users dot sourceforge dot net
  Oliver Ebenhöh

SEE ALSO

Top

New versions of this module can be found on http://steffen-mueller.net or CPAN. The module development takes place on Sourceforge at http://sourceforge.net/projects/math-symbolic/

Math::Symbolic

Math-Symbolic documentation Contained in the Math-Symbolic distribution. 
package Math::Symbolic::Derivative;

use 5.006;
use strict;
use warnings;
no warnings 'recursion';

use Carp;

use Math::Symbolic::ExportConstants qw/:all/;

require Exporter;

our @ISA = qw(Exporter);

our %EXPORT_TAGS = (
    'all' => [
        qw(
          &total_derivative
          &partial_derivative
          )
    ]
);

our @EXPORT_OK = ( @{ $EXPORT_TAGS{'all'} } );

our @EXPORT = qw();

our $VERSION = '0.606';


# lookup-table for derivative rules for various operators.
our %Rules = (
    'each operand'                      => \&_each_operand,
    'product rule'                      => \&_product_rule,
    'quotient rule'                     => \&_quotient_rule,
    'logarithmic chain rule after ln'   => \&_logarithmic_chain_rule_after_ln,
    'logarithmic chain rule'            => \&_logarithmic_chain_rule,
    'derivative commutation'            => \&_derivative_commutation,
    'trigonometric derivatives'         => \&_trigonometric_derivatives,
    'inverse trigonometric derivatives' => \&_inverse_trigonometric_derivatives,
    'inverse atan2'                     => \&_inverse_atan2,
);

# References to derivative subroutines
# Will be assigned a reference after subroutine compilation.
our $Partial_Sub;
our $Total_Sub;

our @Constant_Simplify = (
    # B_SUM
    sub {
        my $tree = shift;
        my ($op1, $op2) = @{$tree->{operands}};
        my ($t1, $t2) = ($op1->term_type(), $op2->term_type());
        if ($t1 == T_CONSTANT) {
            return $op2 if $op1->{value} == 0;
            if ($t2 == T_CONSTANT) {
                return Math::Symbolic::Constant->new($op1->{value} + $op2->{value});
            }
        }
        elsif ($t2 == T_CONSTANT) {
            return $op1 if $op2->{value} == 0;
        }

        return $tree;
    },

    # B_DIFFERENCE
    sub {
        my $tree = shift;
        my ($op1, $op2) = @{$tree->{operands}};
        my ($t1, $t2) = ($op1->term_type(), $op2->term_type());
        if ($t1 == T_CONSTANT) {
            $op2 *= -1, return $op2 if $op1->{value} == 0;
            if ($t2 == T_CONSTANT) {
                return Math::Symbolic::Constant->new($op1->{value} - $op2->{value});
            }
        }
        elsif ($t2 == T_CONSTANT) {
            return $op1 if $op2->{value} == 0;
            $op2->{value} *= -1;
            return Math::Symbolic::Operator->new('+', $op1, $op2);
        }
        return $tree;
    },
    
    # B_PRODUCT
    undef, # implemented inline
    # B_DIVISION
    undef, # not implemented

    # U_MINUS
    sub {
        my $tree = shift;
        my $op = $tree->{operands}[0];
        if ($op->term_type == T_CONSTANT) {
            return Math::Symbolic::Constant->new(-$op->{value});
        }
        return $tree;
    },

    #... not implemented
);


sub _each_operand {
    my ( $tree, $var, $cloned, $d_sub ) = @_;
    foreach ( @{ $tree->{operands} } ) {
        $_ = $d_sub->( $_, $var, 1 );
    }

    my $type = $tree->type();
    my $simplifier = $Constant_Simplify[$type];
    return $simplifier->($tree) if $simplifier;

    return $tree;
}


sub _product_rule {
    my ( $tree, $var, $cloned, $d_sub ) = @_;
    my $ops = $tree->{operands};
    my ($o1, $o2) = @$ops;
    my ($to1, $to2) = ($o1->term_type(), $o2->term_type());

    # one of the terms is a constant, don't derive it
    if ($to1 == T_CONSTANT) {
        return Math::Symbolic::Constant->zero() if $o1->{value} == 0;
        my $deriv = $d_sub->( $o2, $var, 0 );
        return $deriv if $o1->{value} == 0;
        return Math::Symbolic::Constant->new($deriv->{value}*$o1->{value})
          if $deriv->term_type == T_CONSTANT;
    }
    if ($to2 == T_CONSTANT) {
        return Math::Symbolic::Constant->zero() if $o2->{value} == 0;
        my $deriv = $d_sub->( $o1, $var, 0 );
        return $deriv if $o2->{value} == 0;
        return Math::Symbolic::Constant->new($deriv->{value}*$o2->{value})
          if $deriv->term_type == T_CONSTANT;
    }
    
    my $do1 = $d_sub->( $o1, $var, 0 );
    my $do2 = $d_sub->( $o2, $var, 0 );

    my ($tdo1, $tdo2) = ($do1->term_type(), $do2->term_type());

    my ($m1, $m2);
    # check for const*const
    if ($tdo1 == T_CONSTANT) {
        if ($to2 == T_CONSTANT) {
            $m1 = $do1->new($o2->{value} * $do1->{value}); # const
        } elsif ($do1->{value} == 0) {
            $m1 = $do1->zero(); # 0
        } elsif ($do1->{value} == 1) {
            $m1 = $o2;
        } else {
            $m1 = $do1*$o2; # c*tree
        }
    }
    else {
        $m1 = $o2*$do1;
    }

    if ($tdo2 == T_CONSTANT) {
        if ($to1 == T_CONSTANT) {
            $m2 = $do2->new($o1->{value} * $do2->{value}); # const
        } elsif ($do2->{value} == 0) {
            $m2 = $do2->zero(); # 0
        } elsif ($do2->{value} == 1) {
            $m2 = $o1;
        } else {
            $m2 = $do2*$o1; # c*tree
        }
    }
    else {
        $m2 = $o1*$do2;
    }

    # 0's or 2 consts in +
    if ($m1->term_type == T_CONSTANT) {
        return $m2 if $m1->{value} == 0;
        if ($m2->term_type == T_CONSTANT) {
            return $m2->new($m1->{value}*$m2->{value});
        }
    }
    elsif ($m2->term_type == T_CONSTANT) {
        return $m1 if $m2->{value} == 0;
    }

    return Math::Symbolic::Operator->new( '+', $m1, $m2 );
}

sub _quotient_rule {
    my ( $tree, $var, $cloned, $d_sub ) = @_;

    my ($op1, $op2) = @{$tree->{operands}};

    my ($do1, $do2);

    # y = f(x)/c; y' = f'/c
    if ($op2->is_simple_constant()) {
        $do1 = $d_sub->( $op1, $var, 0 );
        my $val = $op2->value();

        if ($val == 0) {
            return $tree->new('/', $do1, $op2->new()); # inf!
        }
        elsif ($val == 1) {
            return $do1; # f/1
        }
        return $tree->new('*', Math::Symbolic::Constant->new(1/$val), $do1);
    }
    # y = c/f(x) => y' = -c*f'(x)/f^2(x)
    elsif ($op1->is_simple_constant()) {
        $do2 = $d_sub->( $op2, $var, 0 );
        my $val = $op1->value();
        
        if ($val == 0) {
            return Math::Symbolic::Constant->zero(); # 0*f'/f
        }

        my $tdo2 = $do2->term_type();
        if ($tdo2 == T_CONSTANT) {
            return $do2->zero() if $do2->{value} == 0; # c*0/f
            return $tree->new(
                '/', $do2->new(-1.*$val*$do2->{value}),
                     $tree->new('^', $op2, 2)
            );
        }
        else {
            return $tree->new(
                '*', Math::Symbolic::Constant->new(-1*$val),
                $tree->new('/', $do2, $tree->new('^', $op2, Math::Symbolic::Constant->new(2)))
            )
        }
    }

    $do1 = $d_sub->( $op1, $var, 0 ) if not $do1;
    $do2 = $d_sub->( $op2, $var, 0 ) if not $do2;

    my $m1  = Math::Symbolic::Operator->new( '*', $do1, $op2 );
    my $m2  = Math::Symbolic::Operator->new( '*', $op1, $do2 );

    # f' = 0
    if ($do1->is_zero()) {
        $m1 = undef;
    }
    # f' = 1
    elsif ($do1->is_one()) {
        $m1 = $op2->new();
    }

    # g' = 0
    if ($do2->is_zero()) {
        $m2 = undef;
    }
    elsif ($do2->is_one()) {
        $m2 = $op1->new();
    }

    my $upper;
    # -g'f / g^2
    if (not defined $m1) {
        # f'=g'=0
        return Math::Symbolic::Constant->zero() if not defined $m2;
        $upper = $tree->new('neg', $m2);
    }
    # f'g / g^2 = f'/g
    elsif (not defined $m2) {
        return $tree->new('/', $do1, $op2);
    }

    my $m3 = $tree->new('^', $op2, Math::Symbolic::Constant->new(2));
    if (not defined $upper) {
      $upper = Math::Symbolic::Operator->new( '-', $m1, $m2 );
    }
    return Math::Symbolic::Operator->new( '/', $upper, $m3 );
}

sub _logarithmic_chain_rule_after_ln {
    my ( $tree, $var, $cloned, $d_sub ) = @_;

    # y(x)=u^v
    # y'(x)=y*(d/dx ln(y))
    # y'(x)=y*(d/dx (v*ln(u)))
    my ($u, $v) = @{$tree->{operands}};

    # This is a special case:
    # y(x)=u^CONST
    # y'(x)=CONST*y* d/dx ln(u)
    # y'(x)=CONST*y* u' / u
    if ($v->term_type() == T_CONSTANT) {

        # y=VAR^CONST
        if ($u->term_type() == T_VARIABLE) {
            my $d = $d_sub->($u, $var, 0);
            my $dtt = $d->term_type();
            if ($dtt == T_CONSTANT) {
                # not our var
                return Math::Symbolic::Constant->zero() if $d->{value} == 0;
                # our var
                return Math::Symbolic::Constant->one() if $v->{value} == 1;
                return $tree->new('*', $v->new(), $u->new()) if $v->{value} == 2;
                return $tree->new('*', $v->new(), $tree->new('^', $u->new(), $v->new($v->{value}-1)));
            }
            # otherwise: signature contains $var
        }
        return Math::Symbolic::Operator->new(
            '*', 
            Math::Symbolic::Operator->new(
                '*', $v->new(), $tree
            ),
            Math::Symbolic::Operator->new(
                '/', $d_sub->($u, $var, 0), $u->new()
            )
        );
    }

    my $e    = Math::Symbolic::Constant->euler();
    my $ln   = Math::Symbolic::Operator->new( 'log', $e, $u );
    my $mul1 = $ln->new( '*', $v, $ln );
    my $dmul = $d_sub->( $mul1, $var, 0 );
    $tree = $ln->new( '*', $tree, $dmul );
    return $tree;
}

sub _logarithmic_chain_rule {
    my ( $tree, $var, $cloned, $d_sub ) = @_;

    #log_a(y(x))=>y'(x)/(ln(a)*y(x))
    my ($a, $y) = @{$tree->{operands}};
    my $dy  = $d_sub->( $y, $var, 0 );

    # This would be y'/y
    if ($a->term_type() == T_CONSTANT and $a->{special} eq 'euler') {
        return Math::Symbolic::Operator->new('/', $dy, $y);
    }
    
    my $e    = Math::Symbolic::Constant->euler();
    my $ln   = Math::Symbolic::Operator->new( 'log', $e, $a );
    my $mul1 = $ln->new( '*', $ln, $y->new() );
    $tree = $ln->new( '/', $dy, $mul1 );
    return $tree;
}

sub _derivative_commutation {
    my ( $tree, $var, $cloned, $d_sub ) = @_;
    $tree->{operands}[0] = $d_sub->( $tree->{operands}[0], $var, 0 );
    return $tree;
}

sub _trigonometric_derivatives {
    my ( $tree, $var, $cloned, $d_sub ) = @_;
    my $op = Math::Symbolic::Operator->new();
    my $d_inner = $d_sub->( $tree->{operands}[0], $var, 0 );
    my $trig;
    my $type = $tree->type();
    if ( $type == U_SINE ) {
        $trig = $op->new( 'cos', $tree->{operands}[0] );
    }
    elsif ( $type == U_COSINE ) {
        $trig = $op->new( 'neg', $op->new( 'sin', $tree->{operands}[0] ) );
    }
    elsif ( $type == U_SINE_H ) {
        $trig = $op->new( 'cosh', $tree->{operands}[0] );
    }
    elsif ( $type == U_COSINE_H ) {
        $trig = $op->new( 'sinh', $tree->{operands}[0] );
    }
    elsif ( $type == U_TANGENT or $type == U_COTANGENT ) {
        $trig = $op->new(
            '/',
            Math::Symbolic::Constant->one(),
            $op->new(
                '^',
                $op->new( 'cos', $tree->op1() ),
                Math::Symbolic::Constant->new(2)
            )
        );
        $trig = $op->new( 'neg', $trig ) if $type == U_COTANGENT;
    }
    else {
        die "Trigonometric derivative applied to invalid operator.";
    }
    if ($d_inner->term_type() == T_CONSTANT) {
        my $spec = $d_inner->special();
        if ($spec eq 'zero') {
            return $d_inner;
        }
        elsif ($spec eq 'one') {
            return $trig;
        }
    }
    return $op->new( '*', $d_inner, $trig );
}

sub _inverse_trigonometric_derivatives {
    my ( $tree, $var, $cloned, $d_sub ) = @_;
    my $op = Math::Symbolic::Operator->new();
    my $d_inner = $d_sub->( $tree->{operands}[0], $var, 0 );
    my $trig;
    my $type = $tree->type();
    if ( $type == U_ARCSINE or $type == U_ARCCOSINE ) {
        my $one = $type == U_ARCSINE
            ? Math::Symbolic::Constant->one()
            : Math::Symbolic::Constant->new(-1);
        $trig = $op->new( '/', $one,
            $op->new( '-', $one->new(1), $op->new( '^', $tree->op1(), $one->new(2) ) )
        );
    }
    elsif ($type == U_ARCTANGENT
        or $type == U_ARCCOTANGENT )
    {
        my $one = $type == U_ARCTANGENT
            ? Math::Symbolic::Constant->one()
            : Math::Symbolic::Constant->new(-1);
        $trig = $op->new( '/', $one,
            $op->new( '+', $one->new(1), $op->new( '^', $tree->op1(), $one->new(2) ) )
        );
    }
    elsif ($type == U_AREASINE_H
        or $type == U_AREACOSINE_H )
    {
        my $one = Math::Symbolic::Constant->one();
        $trig = $op->new(
            '/', $one,
            $op->new(
                '^',
                $op->new(
                    ( $tree->type() == U_AREASINE_H ? '+' : '-' ),
                    $op->new( '^', $tree->op1(), $one->new(2) ),
                    $one
                ),
                $one->new(0.5)
            )
        );
    }
    else {
        die "Inverse trig. derivative applied to invalid operator.";
    }

    if ($d_inner->term_type() == T_CONSTANT) {
        my $spec = $d_inner->special();
        if ($spec eq 'zero') {
            return $d_inner;
        }
        elsif ($spec eq 'one') {
            return $trig;
        }
    }
    return $op->new( '*', $d_inner, $trig );
}

sub _inverse_atan2 {
    my ( $tree, $var, $cloned, $d_sub ) = @_;
    # d/df atan(y/x) = x^2/(x^2+y^2) * (d/df y/x)
    my ($op1, $op2) = @{$tree->{operands}}; 

    my $inner = $d_sub->( $op1->new()/$op2->new(), $var, 0 );
    # templates
    my $two = Math::Symbolic::Constant->new(2);
    my $op  = Math::Symbolic::Operator->new('+', $two, $two);

    my $result = $op->new('*',
      $op->new('/',
        $op->new('^', $op2->new(), $two->new()), 
        $op->new(
          '+', $op->new('^', $op2->new(), $two->new()),
          $op->new('^', $op1->new(), $two->new())
        )
      ),
      $inner
    );
    return $result;
}


sub partial_derivative {
    my $tree = shift;
    my $var  = shift;
    defined $var or die "Cannot derive using undefined variable.";
    if ( ref($var) eq '' ) {
        $var = Math::Symbolic::parse_from_string($var);
        croak "2nd argument to partial_derivative must be variable."
          if ( ref($var) ne 'Math::Symbolic::Variable' );
    }
    else {
        croak "2nd argument to partial_derivative must be variable."
          if ( ref($var) ne 'Math::Symbolic::Variable' );
    }

    my $cloned = shift;

    if ( not $cloned ) {
        $tree   = $tree->new();
        $cloned = 1;
    }

    if ( $tree->term_type() == T_OPERATOR ) {
        my $rulename =
          $Math::Symbolic::Operator::Op_Types[ $tree->type() ]->{derive};
        my $subref = $Rules{$rulename};

        die "Cannot derive using rule '$rulename'."
          unless defined $subref;
        $tree = $subref->( $tree, $var, $cloned, $Partial_Sub );
    }
    elsif ( $tree->term_type() == T_CONSTANT ) {
        $tree = Math::Symbolic::Constant->zero();
    }
    elsif ( $tree->term_type() == T_VARIABLE ) {
        if ( $tree->name() eq $var->name() ) {
            $tree = Math::Symbolic::Constant->one;
        }
        else {
            $tree = Math::Symbolic::Constant->zero;
        }
    }
    else {
        die "Cannot apply partial derivative to anything but a tree.";
    }

    return $tree;
}


sub total_derivative {
    my $tree = shift;
    my $var  = shift;
    defined $var or die "Cannot derive using undefined variable.";
    if ( ref($var) eq '' ) {
        $var = Math::Symbolic::parse_from_string($var);
        croak "Second argument to total_derivative must be variable."
          if ( ref($var) ne 'Math::Symbolic::Variable' );
    }
    else {
        croak "Second argument to total_derivative must be variable."
          if ( ref($var) ne 'Math::Symbolic::Variable' );
    }

    my $cloned = shift;

    if ( not $cloned ) {
        $tree   = $tree->new();
        $cloned = 1;
    }

    if ( $tree->term_type() == T_OPERATOR ) {
        my $var_name = $var->name();
        my @tree_sig = $tree->signature();
        if ( ( grep { $_ eq $var_name } @tree_sig ) > 0 ) {
            my $rulename =
              $Math::Symbolic::Operator::Op_Types[ $tree->type() ]->{derive};
            my $subref = $Rules{$rulename};

            die "Cannot derive using rule '$rulename'."
              unless defined $subref;
            $tree = $subref->( $tree, $var, $cloned, $Total_Sub );
        }
        else {
            $tree = Math::Symbolic::Constant->zero();
        }
    }
    elsif ( $tree->term_type() == T_CONSTANT ) {
        $tree = Math::Symbolic::Constant->zero();
    }
    elsif ( $tree->term_type() == T_VARIABLE ) {
        my $name     = $tree->name();
        my $var_name = $var->name();

        if ( $name eq $var_name ) {
            $tree = Math::Symbolic::Constant->one;
        }
        else {
            my @tree_sig = $tree->signature();
            my $is_dependent;
            foreach my $ident (@tree_sig) {
                if ( $ident eq $var_name ) {
                    $is_dependent = 1;
                    last;
                }
            }
            if ( $is_dependent ) {
                $tree =
                  Math::Symbolic::Operator->new( 'total_derivative', $tree,
                    $var );
            }
            else {
                $tree = Math::Symbolic::Constant->zero;
            }
        }
    }
    else {
        die "Cannot apply total derivative to anything but a tree.";
    }

    return $tree;
}

# Class data again.
$Partial_Sub = \&partial_derivative;
$Total_Sub   = \&total_derivative;

1;
__END__