login  home  contents  what's new  discussion  bug reports help  links  subscribe  changes  refresh  edit

# Edit detail for SymbolicDifferentiation revision 13 of 13

 1 2 3 4 5 6 7 8 9 10 11 12 13 Editor: hemmecke Time: 2014/02/27 12:34:27 GMT+0 Note:

removed:
-    dprint a ==> rhx:="rhx" -- print a

removed:
-        dprint(f "plusD: " + f l)

removed:
-        dprint(f "timesD: " + f l)

removed:
-        dprint(f "timesD len: " + f len)

removed:
-            dprint(f "dTimes i: " + f i)

removed:
-                dprint(f "timesD j: " + f j)
-                dprint(f "timesD z: " + f z)

removed:
-        dprint(f "powerD: " + f l)

removed:
-        dprint(f "powerXD: " + f l)

removed:
-        dprint(f "nthRootD: " + f l)


## Symbolic Differentiation Tutorial

In the following, we like to program a simple machinery that allows us to differentiate univariate functions.

To understand the programs here, you should first have read and understood ProgrammingSPAD.

The aim of this page is rather to give a short account on how one can quickly achieve quite some sophisticated program with a few lines of SPAD code. The code deliberately does only build on a very small subset of the FriCAS library. In fact, we only use the category SetCategory and the domains Integer, Fraction, Complex, String, Symbol, List, Product, Kernel, OutputForm, Expression. Although the Expression domain already comes with differentiation, we are not going to use it. We could have done without Expression at the cost of implementing our own expression domain rule for identifying e1+2 with e2+e1, etc. However, since the purpose of this page is to show programming in SPAD, rather than entering the complicated area of low-level expression manipulation, we decided against it.

The first file mex.spad was actually used as the basis for a little project to teach some children (age 16-18) how Symbolic Differentiation basically works. They were supposed to "invent" the dex.spad file.

The file mex.spad is not really dealing with symbolic differentiation. It rather deals with very basic things like defining what a monoid, a ring, or a field is.

The most interesting thing is that it defines expressions that can be decomposed into a top-level symbol and its corresponding arguments (which might be complicated expressions themselves).

With the help of this MExpression domain, it is realatively easy to program symbolic differentiation by simply looking at the top-level symbol, and applying the corresponding rule for it and differentiating the arguments according to the chain rule.

fricas
(1) -> <spad>
-------------------------------------------------------------------
--
-- MEX
-- Copyright (C) 2012,  Ralf Hemmecke <ralf@hemmecke.de>
--
-------------------------------------------------------------------
-- This program is free software: you can redistribute it and/or modify
-- the Free Software Foundation, either version 3 of the License, or
-- (at your option) any later version.
--
-- This program is distributed in the hope that it will be useful,
-- but WITHOUT ANY WARRANTY; without even the implied warranty of
-- MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
-- GNU General Public License for more details.
--
-- You should have received a copy of the GNU General Public License
-- along with this program.  If not, see <http://www.gnu.org/licenses/>.
-------------------------------------------------------------------
-- This file contains basic definitions for a simple expression domain.
-------------------------------------------------------------------
-- These two macros are necessary to distinguish between Rep and %.
rep x ==> (x@%) pretend Rep
per x ==> (x@Rep) pretend %
Z ==> Integer
Q ==> Fraction Z
C ==> Complex Q
X ==> Expression C
D ==> Product(String, List %)
-------------------------------------------------------------------
fricas
)abbrev category MABMON MAbelianAdditiveMonoid
++ MAbelianAdditiveMonoid is the category of abelian monoids that
++ are written in an additive way.
0: %
++ 0 is the neutral element w.r.t. +.
_+: (%, %) -> %
++ + is a commutative and associative operation.
-------------------------------------------------------------------
fricas
)abbrev category MABAGRP MAbelianAdditiveGroup
++ MAbelianAdditiveGroup is the category of abelian groups that
++ are written in an additive way.
_-: % -> %
++ x + (-x) == 0
_-: (%, %) -> %
++ x - y = z if and only if x = y + z.
((x: %) - (y: %)): % == x + (- y)
-------------------------------------------------------------------
fricas
)abbrev category MMMON MMultiplicativeMonoid
++ MMultiplicativeMonoid is the category of monoids that
++ are written in a multiplicative way.
MMultiplicativeMonoid: Category == SetCategory with
1: %
++ 1 is the neutral element w.r.t. *.
_*: (%, %) -> %
++ * is an associative operation.
-------------------------------------------------------------------
fricas
)abbrev category MABMMON MAbelianMultiplicativeMonoid
++ MAbelianMultiplicativeMonoid is the category of abelian monoids that
++ are written in a multiplicative way.
MAbelianMultiplicativeMonoid: Category == MMultiplicativeMonoid
-------------------------------------------------------------------
fricas
)abbrev category MRING MRing
++ MRing is the category of rings with unity.
-------------------------------------------------------------------
fricas
)abbrev category MCRING MCommutativeRing
++ MCommutativeRing is the category of commutative rings with unity.
MCommutativeRing: Category == Join(MRing, MAbelianMultiplicativeMonoid)
-------------------------------------------------------------------
fricas
)abbrev category MINTDOM MIntegralDomain
++ MIntegralDomain is the category of commutative rings with unity
++ that have no zero-divisors.
MIntegralDomain: Category == MCommutativeRing
-------------------------------------------------------------------
fricas
)abbrev category MFIELD MField
++ MField is the category of fields.
MField: Category == MIntegralDomain with
inv: % -> %
++ x * inv(x) = 1
_/: (%, %) -> %
++ Only defined if the second argument is non-zero.
++ x/y = z if and only if x = z*y
_^: (%, Integer) -> %
++ x^n is an abbreviation for x*x*...*x (n-fold product of x) if n>0.
++ x^0 := 1 for any x (even for x=0),
++ x^n := 1/(x^(-n)) if n<0.
inv(x: %): % == 1/x
-------------------------------------------------------------------
fricas
)abbrev category MEXCAT MExpressionCategory
++ MExpressionCategory is the category of formal expressions domains
++ that are closed under field operations and taking sin, cos, exp,
++ log, and sqrt.
MExpressionCategory: Category == MField with
coerce: % -> X
++ coerce(z) transforms z into an equivalent value in
++ Expression(Complex Fraction Integer).
coerce: Z -> %
++ coerce(z) transforms the integer z into an element of this type.
coerce: Q -> %
++ coerce(z) transforms the rational number z into an element
++ of this type.
coerce: Complex Z -> %
++ coerce(z) transforms the complex number z into an element
++ of this type.
coerce: Complex Q -> %
++ coerce(z) transforms the complex number z into an element
++ of this type.
coerce: Symbol -> %
++ coerce(z) returns z as a variable of this type.
variable?: % -> Boolean
++ variable?(z) returns true if z is a variable and false otherwise.
complex?: % -> Boolean
++ complex?(z) returns true if z can be represented as a
++ complex number and false otherwise
_=: (%, %) -> Boolean
++ y=z tests whether two expressions are (after mild simplification)
++ syntactically equal.
decompose: % -> Product(String, List %)
++ decompose(z) return ["complex", [z]] if complex?(z) is true.
++ It returns ["variable", [z]] if variable?(z) is true.
++ In all other cases z is an expression of the form foo(a1,...,an)
++ and ["foo", [a1, ..., an]] will be returned.
++ Note that for infix expression like a1+a2+...+an, this functions
++ returns ["+",[a1,...,an]].
_^: (%, %) -> %
++ y^z is a formal exponentiation expression.
sqrt: % -> %
++ sqrt(z) is an expresssion that denotes square root of z.
sin: % -> %
++ sin(z) is an expression that denotes the sine of z.
cos: % -> %
++ cos(z) is an expression that denotes the cosine of z.
exp: % -> %
++ exp(z) is an expression that denotes the e^z where e is
++ Euler's constant.
log: % -> %
++ log(z) is an expression that denotes the natural logarithm of z.
tan: % -> %
++ tan(z) is an expression that denotes the tangent of z.
cot: % -> %
++ cot(z) is an expression that denotes the cotangent of z.
asin: % -> %
++ asin(z) is an expression that denotes the arc-sine of z.
acos: % -> %
++ acos(z) is an expression that denotes the arc-cosine of z.
atan: % -> %
++ atan(z) is an expression that denotes the arc-tangent of z.
acot: % -> %
++ acot(z) is an expression that denotes the arc-cotangent of z.
-------------------------------------------------------------------
fricas
)abbrev domain MEX MExpression
Rep ==> X
-- exported functions
coerce(a: %): X == a pretend X
coerce(a: %): OutputForm == coerce rep a
coerce(z: Z): % == per(z::X)
coerce(q: Q): % == per(q::X)
coerce(c: C): % == per(c::X)
coerce(c: Complex Z): % == complex(real(c)::Q, imag(c)::Q)$C::% coerce(s: Symbol): % == per(s::X) variable?(x: %): Boolean == empty? second decompose x complex?(a: %): Boolean == number? rep a 0: % == (0$Z)::%
1: % == (1$Z)::% ((a: %) = (b: %)): Boolean == rep a = rep b _-(a: %): % == per(- rep a) ((a: %) + (b: %)): % == per(rep a + rep b) ((a: %) * (b: %)): % == per(rep a * rep b) ((a: %) / (b: %)): % == per(rep a / rep b) ((a: %) ^ (z: Z)): % == per(rep a ^ z) ((a: %) ^ (b: %)): % == per(rep a ^ rep b) sqrt(a: %): % == per sqrt rep a sin(a: %): % == per sin rep a cos(a: %): % == per cos rep a exp(a: %): % == per exp rep a log(a: %): % == per log rep a tan(a: %): % == per tan rep a cot(a: %): % == per cot rep a asin(a: %): % == per asin rep a acos(a: %): % == per acos rep a atan(a: %): % == per atan rep a acot(a: %): % == per acot rep a decompose(a: %): D == complex? a => ["complex", [a]]$D
ra := rep a
il := isPlus ra
il case List(X) => ["+", (il::List(X)) pretend List(%)]$D il := isTimes ra il case List(X) => ["*", (il::List(X)) pretend List(%)]$D
ip := isPower ra
ip case Record(val: X, exponent: Z) and ip.exponent ~= 1 =>
["^", [ip.val, (ip.exponent)::X] pretend List(%)]$D ks: List(Kernel X) := kernels ra #(ks)~=1 => print("decompose: more than one kernel detected"::OutputForm) print(ra::OutputForm) print(ks::OutputForm) error "decompose error" k: Kernel(X) := first ks [string name k, [per z for z in argument k]]$D</spad>
fricas
Compiling FriCAS source code from file
using old system compiler.
------------------------------------------------------------------------
compiling into NRLIB MABMON
Time: 0.01 SEC.
finalizing NRLIB MABMON
--------constructor---------
--------((Zero) (%) constant)---------
--------(+ (% % %))---------
; compiling file "/var/aw/var/LatexWiki/MABMON.NRLIB/MABMON.lsp" (written 21 FEB 2023 12:47:03 PM):
; /var/aw/var/LatexWiki/MABMON.NRLIB/MABMON.fasl written
; compilation finished in 0:00:00.006
------------------------------------------------------------------------
MAbelianAdditiveMonoid is now explicitly exposed in frame initial
/var/aw/var/LatexWiki/MABMON.NRLIB/MABMON
------------------------------------------------------------------------
compiling into NRLIB MABAGRP
Time: 0 SEC.
------------------------------------------------------------------------
compiling into NRLIB MABAGRP-
compiling exported - : (S,S) -> S
Time: 0 SEC.
(time taken in buildFunctor:  0)
Time: 0 SEC.
Time: 0 seconds
finalizing NRLIB MABAGRP-
--------constructor---------
--------(- (% %))---------
--------(- (% % %))---------
; compiling file "/var/aw/var/LatexWiki/MABAGRP-.NRLIB/MABAGRP-.lsp" (written 21 FEB 2023 12:47:03 PM):
; /var/aw/var/LatexWiki/MABAGRP-.NRLIB/MABAGRP-.fasl written
; compilation finished in 0:00:00.011
------------------------------------------------------------------------
MAbelianAdditiveGroup& is now explicitly exposed in frame initial
/var/aw/var/LatexWiki/MABAGRP-.NRLIB/MABAGRP-
finalizing NRLIB MABAGRP
--------constructor---------
--------(- (% %))---------
--------(- (% % %))---------
; compiling file "/var/aw/var/LatexWiki/MABAGRP.NRLIB/MABAGRP.lsp" (written 21 FEB 2023 12:47:03 PM):
; /var/aw/var/LatexWiki/MABAGRP.NRLIB/MABAGRP.fasl written
; compilation finished in 0:00:00.003
------------------------------------------------------------------------
MAbelianAdditiveGroup is now explicitly exposed in frame initial
/var/aw/var/LatexWiki/MABAGRP.NRLIB/MABAGRP
MMMON abbreviates category MMultiplicativeMonoid
------------------------------------------------------------------------
initializing NRLIB MMMON for MMultiplicativeMonoid
compiling into NRLIB MMMON
;;;     ***       |MMultiplicativeMonoid| REDEFINED
Time: 0 SEC.
finalizing NRLIB MMMON
Processing MMultiplicativeMonoid for Browser database:
--------constructor---------
--------((One) (%) constant)---------
--->-->MMultiplicativeMonoid(((One) (%) constant)): Improper first word in comments:
--------(* (% % %))---------
; compiling file "/var/aw/var/LatexWiki/MMMON.NRLIB/MMMON.lsp" (written 21 FEB 2023 12:47:03 PM):
; /var/aw/var/LatexWiki/MMMON.NRLIB/MMMON.fasl written
; compilation finished in 0:00:00.002
------------------------------------------------------------------------
MMultiplicativeMonoid is now explicitly exposed in frame initial
MMultiplicativeMonoid will be automatically loaded when needed from
/var/aw/var/LatexWiki/MMMON.NRLIB/MMMON
MABMMON abbreviates category MAbelianMultiplicativeMonoid
------------------------------------------------------------------------
initializing NRLIB MABMMON for MAbelianMultiplicativeMonoid
compiling into NRLIB MABMMON
;;;     ***       |MAbelianMultiplicativeMonoid| REDEFINED
Time: 0 SEC.
finalizing NRLIB MABMMON
Processing MAbelianMultiplicativeMonoid for Browser database:
--------constructor---------
; compiling file "/var/aw/var/LatexWiki/MABMMON.NRLIB/MABMMON.lsp" (written 21 FEB 2023 12:47:03 PM):
; /var/aw/var/LatexWiki/MABMMON.NRLIB/MABMMON.fasl written
; compilation finished in 0:00:00.003
------------------------------------------------------------------------
MAbelianMultiplicativeMonoid is now explicitly exposed in frame
initial
MAbelianMultiplicativeMonoid will be automatically loaded when
needed from /var/aw/var/LatexWiki/MABMMON.NRLIB/MABMMON
MRING abbreviates category MRing
------------------------------------------------------------------------
initializing NRLIB MRING for MRing
compiling into NRLIB MRING
;;;     ***       |MRing| REDEFINED
Time: 0 SEC.
finalizing NRLIB MRING
Processing MRing for Browser database:
--------constructor---------
; compiling file "/var/aw/var/LatexWiki/MRING.NRLIB/MRING.lsp" (written 21 FEB 2023 12:47:03 PM):
; /var/aw/var/LatexWiki/MRING.NRLIB/MRING.fasl written
; compilation finished in 0:00:00.002
------------------------------------------------------------------------
MRing is now explicitly exposed in frame initial
MRing will be automatically loaded when needed from
/var/aw/var/LatexWiki/MRING.NRLIB/MRING
MCRING abbreviates category MCommutativeRing
------------------------------------------------------------------------
initializing NRLIB MCRING for MCommutativeRing
compiling into NRLIB MCRING
;;;     ***       |MCommutativeRing| REDEFINED
Time: 0 SEC.
finalizing NRLIB MCRING
Processing MCommutativeRing for Browser database:
--------constructor---------
; compiling file "/var/aw/var/LatexWiki/MCRING.NRLIB/MCRING.lsp" (written 21 FEB 2023 12:47:03 PM):
; /var/aw/var/LatexWiki/MCRING.NRLIB/MCRING.fasl written
; compilation finished in 0:00:00.002
------------------------------------------------------------------------
MCommutativeRing is now explicitly exposed in frame initial
MCommutativeRing will be automatically loaded when needed from
/var/aw/var/LatexWiki/MCRING.NRLIB/MCRING
MINTDOM abbreviates category MIntegralDomain
------------------------------------------------------------------------
initializing NRLIB MINTDOM for MIntegralDomain
compiling into NRLIB MINTDOM
;;;     ***       |MIntegralDomain| REDEFINED
Time: 0 SEC.
finalizing NRLIB MINTDOM
Processing MIntegralDomain for Browser database:
--------constructor---------
; compiling file "/var/aw/var/LatexWiki/MINTDOM.NRLIB/MINTDOM.lsp" (written 21 FEB 2023 12:47:03 PM):
; /var/aw/var/LatexWiki/MINTDOM.NRLIB/MINTDOM.fasl written
; compilation finished in 0:00:00.003
------------------------------------------------------------------------
MIntegralDomain is now explicitly exposed in frame initial
MIntegralDomain will be automatically loaded when needed from
/var/aw/var/LatexWiki/MINTDOM.NRLIB/MINTDOM
MFIELD abbreviates category MField
------------------------------------------------------------------------
initializing NRLIB MFIELD for MField
compiling into NRLIB MFIELD
;;;     ***       |MField| REDEFINED
Time: 0 SEC.
MFIELD- abbreviates domain MField&
------------------------------------------------------------------------
initializing NRLIB MFIELD- for MField&
compiling into NRLIB MFIELD-
compiling exported inv : S -> S
Time: 0.01 SEC.
(time taken in buildFunctor:  0)
;;;     ***       |MField&| REDEFINED
Time: 0 SEC.
Cumulative Statistics for Constructor MField&
Time: 0.01 seconds
finalizing NRLIB MFIELD-
Processing MField& for Browser database:
--------constructor---------
--------(inv (% %))---------
--->-->MField&((inv (% %))): Improper initial operator in comments: *
--------(/ (% % %))---------
--->-->MField&((/ (% % %))): Improper first word in comments: Only
"Only defined if the second argument is non-zero. x/y = \\spad{z} if and only if \\spad{x} = z*y"
--------(^ (% % (Integer)))---------
; compiling file "/var/aw/var/LatexWiki/MFIELD-.NRLIB/MFIELD-.lsp" (written 21 FEB 2023 12:47:03 PM):
; /var/aw/var/LatexWiki/MFIELD-.NRLIB/MFIELD-.fasl written
; compilation finished in 0:00:00.008
------------------------------------------------------------------------
MField& is now explicitly exposed in frame initial
MField& will be automatically loaded when needed from
/var/aw/var/LatexWiki/MFIELD-.NRLIB/MFIELD-
finalizing NRLIB MFIELD
Processing MField for Browser database:
--------constructor---------
--------(inv (% %))---------
--->-->MField((inv (% %))): Improper initial operator in comments: *
--------(/ (% % %))---------
--->-->MField((/ (% % %))): Improper first word in comments: Only
"Only defined if the second argument is non-zero. x/y = \\spad{z} if and only if \\spad{x} = z*y"
--------(^ (% % (Integer)))---------
; compiling file "/var/aw/var/LatexWiki/MFIELD.NRLIB/MFIELD.lsp" (written 21 FEB 2023 12:47:03 PM):
; /var/aw/var/LatexWiki/MFIELD.NRLIB/MFIELD.fasl written
; compilation finished in 0:00:00.002
------------------------------------------------------------------------
MField is now explicitly exposed in frame initial
MField will be automatically loaded when needed from
/var/aw/var/LatexWiki/MFIELD.NRLIB/MFIELD
MEXCAT abbreviates category MExpressionCategory
------------------------------------------------------------------------
initializing NRLIB MEXCAT for MExpressionCategory
compiling into NRLIB MEXCAT
;;;     ***       |MExpressionCategory| REDEFINED
Time: 0 SEC.
finalizing NRLIB MEXCAT
Processing MExpressionCategory for Browser database:
--------constructor---------
--------(coerce ((Expression (Complex (Fraction (Integer)))) %))---------
--------(coerce (% (Integer)))---------
--------(coerce (% (Fraction (Integer))))---------
--------(coerce (% (Complex (Integer))))---------
--------(coerce (% (Complex (Fraction (Integer)))))---------
--------(coerce (% (Symbol)))---------
--------(variable? ((Boolean) %))---------
--------(complex? ((Boolean) %))---------
--------(= ((Boolean) % %))---------
--------(decompose ((Product (String) (List %)) %))---------
--------(^ (% % %))---------
--------(sqrt (% %))---------
--------(sin (% %))---------
--------(cos (% %))---------
--------(exp (% %))---------
--------(log (% %))---------
--------(tan (% %))---------
--------(cot (% %))---------
--------(asin (% %))---------
--------(acos (% %))---------
--------(atan (% %))---------
--------(acot (% %))---------
; compiling file "/var/aw/var/LatexWiki/MEXCAT.NRLIB/MEXCAT.lsp" (written 21 FEB 2023 12:47:03 PM):
; /var/aw/var/LatexWiki/MEXCAT.NRLIB/MEXCAT.fasl written
; compilation finished in 0:00:00.003
------------------------------------------------------------------------
MExpressionCategory is now explicitly exposed in frame initial
MExpressionCategory will be automatically loaded when needed from
/var/aw/var/LatexWiki/MEXCAT.NRLIB/MEXCAT
MEX abbreviates domain MExpression
------------------------------------------------------------------------
initializing NRLIB MEX for MExpression
compiling into NRLIB MEX
processing macro definition Rep ==> Expression Complex Fraction Integer
compiling exported coerce : $-> Expression Complex Fraction Integer MEX;coerce;$E;1 is replaced by a
Time: 0.06 SEC.
compiling exported coerce : $-> OutputForm Time: 0 SEC. compiling exported coerce : Integer ->$
Time: 0 SEC.
compiling exported coerce : Fraction Integer -> $Time: 0 SEC. compiling exported coerce : Complex Fraction Integer ->$
Time: 0.01 SEC.
compiling exported coerce : Complex Integer -> $Time: 0 SEC. compiling exported coerce : Symbol ->$
Time: 0 SEC.
compiling exported variable? : $-> Boolean Time: 0.02 SEC. compiling exported complex? :$ -> Boolean
Time: 0.01 SEC.
compiling exported Zero : () -> $Time: 0 SEC. compiling exported One : () ->$
Time: 0 SEC.
compiling exported = : ($,$) -> Boolean
Time: 0 SEC.
compiling exported - : $->$
Time: 0 SEC.
compiling exported + : ($,$) -> $Time: 0 SEC. compiling exported * : ($,$) ->$
Time: 0 SEC.
compiling exported / : ($,$) -> $Time: 0.01 SEC. compiling exported ^ : ($,Integer) -> $Time: 0 SEC. compiling exported ^ : ($,$) ->$
Time: 0.01 SEC.
compiling exported sqrt : $->$
Time: 0 SEC.
compiling exported sin : $->$
Time: 0.01 SEC.
compiling exported cos : $->$
Time: 0 SEC.
compiling exported exp : $->$
Time: 0.01 SEC.
compiling exported log : $->$
Time: 0 SEC.
compiling exported tan : $->$
Time: 0 SEC.
compiling exported cot : $->$
Time: 0 SEC.
compiling exported asin : $->$
Time: 0 SEC.
compiling exported acos : $->$
Time: 0.01 SEC.
compiling exported atan : $->$
Time: 0 SEC.
compiling exported acot : $->$
Time: 0 SEC.
compiling exported decompose : $-> Product(String,List$)
Time: 0.06 SEC.
(time taken in buildFunctor:  0)
;;;     ***       |MExpression| REDEFINED
;;;     ***       |MExpression| REDEFINED
Time: 0.01 SEC.
Warnings:
[1] variable?:   has no value
[2] decompose:  exponent has no value
[3] decompose:  val has no value
Cumulative Statistics for Constructor MExpression
Time: 0.22 seconds
finalizing NRLIB MEX
Processing MExpression for Browser database:
--->-->MExpression(): Missing Description
; compiling file "/var/aw/var/LatexWiki/MEX.NRLIB/MEX.lsp" (written 21 FEB 2023 12:47:03 PM):
; /var/aw/var/LatexWiki/MEX.NRLIB/MEX.fasl written
; compilation finished in 0:00:00.075
------------------------------------------------------------------------
MExpression is now explicitly exposed in frame initial
MExpression will be automatically loaded when needed from
/var/aw/var/LatexWiki/MEX.NRLIB/MEX

FILE ==> "dex.spad"
-------------------------------------------------------------------
--
-- DEX
-- Copyright (C) 2012,  Ralf Hemmecke <ralf@hemmecke.de>
--
-------------------------------------------------------------------
-- This program is free software: you can redistribute it and/or modify
-- the Free Software Foundation, either version 3 of the License, or
-- (at your option) any later version.
--
-- This program is distributed in the hope that it will be useful,
-- but WITHOUT ANY WARRANTY; without even the implied warranty of
-- MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
-- GNU General Public License for more details.
--
-- You should have received a copy of the GNU General Public License
-- along with this program.  If not, see <http://www.gnu.org/licenses/>.
-------------------------------------------------------------------
-- This file extends a basic expression domain by differentation.
-------------------------------------------------------------------
)abbrev domain DEX DExpression
D ==> Product(String, List %)
DExpression: MExpressionCategory with
differentiate: (%, %) -> %
f a ==> a::OutputForm
-- auxiliary functions
plusD(l: List %, x: %): % ==
r: % := 0
for z in l repeat r := r +  differentiate(z, x)
r
timesD(l: List %, x: %): % ==
s: % := 0
len := #(l)
for i in 1..len repeat
p: % := 1
for z in l for j in 1..len repeat
p := p*(if i=j then differentiate(z, x) else z)
s := s + p
s
powerD(l: List %, x: %): % ==
#(l) ~= 2 => error "powerD: not exactly 2 arguments"
z: % := first l
e: % := first rest l
not complex? e => error "powerD: second argument is not complex"
e*(z^(e-1))*differentiate(z, x)
powerXD(l: List %, x: %): % ==
#(l) ~= 2 => error "powerXD: not exactly 2 arguments"
z: % := first l
e: % := first rest l
complex? z => log z * z^e * differentiate(e, x)
differentiate(exp(log z * e), x)
chainD(l: List %, x: %, d: % -> %): % ==
#(l) ~= 1 => error "chainD: not exactly one argument"
z := first l
d(z)*differentiate(z, x)
nthRootD(l: List %, x: %): % ==
#(l) ~= 2 => error "nthRootD: not exactly two arguments"
z := first l
e := inv first rest l
not complex? e => error "nthRootD: second argument is not complex"
e*(z^(e-1))*differentiate(z, x)
-- exported functions
differentiate(z: %, x: %): % ==
not variable? x => error "second argument must be a variable"
complex? z => 0
z=x => 1
variable? z => 0
d: D := decompose z
o: String := first d
a: List % := second d
o = "+" => plusD(a, x)
o = "*" => timesD(a, x)
o = "^" => powerD(a, x)
o = "%power" => powerXD(a, x)
o = "nthRoot" => nthRootD(a, x)
o = "sin" => chainD(a, x, cos)
o = "cos" => chainD(a, x, t +-> - sin t)
o = "exp" => chainD(a, x, exp)
o = "log" => chainD(a, x, t +-> 1/t)
o = "tan" => chainD(a, x, t +-> tan(t)^2+1)
o = "cot" => chainD(a, x, t +-> -cot(t)^2-1)
o = "asin" => chainD(a, x, t +-> 1/sqrt(1-t^2))
o = "acos" => chainD(a, x, t +-> -1/sqrt(1-t^2))
o = "atan" => chainD(a, x, t +-> 1/(1+t^2))
o = "acot" => chainD(a, x, t +-> -1/(1+t^2))
print("differentiate expression: "::OutputForm)
print(z::OutputForm)
print(((z=x)@Boolean)::OutputForm)
print(d::OutputForm)
print(o::OutputForm)
print(a::OutputForm)
error "Cannot differentiate"
   Compiling FriCAS source code from file
using old system compiler.
DEX abbreviates domain DExpression
------------------------------------------------------------------------
initializing NRLIB DEX for DExpression
compiling into NRLIB DEX
processing macro definition f a ==> ::(a,OutputForm)
compiling local plusD : (List $,$) -> $Time: 0 SEC. compiling local timesD : (List$,$) ->$
Time: 0.01 SEC.
compiling local powerD : (List $,$) -> $Time: 0 SEC. compiling local powerXD : (List$,$) ->$
Time: 0.01 SEC.
compiling local chainD : (List $,$,$->$) -> $Time: 0 SEC. compiling local nthRootD : (List$,$) ->$
Time: 0 SEC.
compiling exported differentiate : ($,$) -> $Time: 0.01 SEC. (time taken in buildFunctor: 0) ;;; *** |DExpression| REDEFINED ;;; *** |DExpression| REDEFINED Time: 0.03 SEC. Cumulative Statistics for Constructor DExpression Time: 0.06 seconds finalizing NRLIB DEX Processing DExpression for Browser database: --->-->DExpression(constructor): Not documented!!!! --->-->DExpression((differentiate (% % %))): Not documented!!!! --->-->DExpression(): Missing Description ; compiling file "/var/aw/var/LatexWiki/DEX.NRLIB/DEX.lsp" (written 21 FEB 2023 12:47:03 PM): ; /var/aw/var/LatexWiki/DEX.NRLIB/DEX.fasl written ; compilation finished in 0:00:00.162 ------------------------------------------------------------------------ DExpression is now explicitly exposed in frame initial DExpression will be automatically loaded when needed from /var/aw/var/LatexWiki/DEX.NRLIB/DEX After having defined an expression domain and having equipped it with a differentiate function, we can now use it. Let's start with some abbreviations. fricas DEX ==> DExpression Type: Void fricas x := 'x :: DEX  (1) Type: DExpression? fricas d(e) ==> differentiate(e, x)$DEX
Type: Void

Now we can differentiate any polynomial formed with the x from above.

fricas
d x
 (2)
Type: DExpression?
fricas
d(x^2)
 (3)
Type: DExpression?
fricas
d(5*(x^2-2)^2)
 (4)
Type: DExpression?

Of course, we can also differentiate all the functions that are included in the code of differentiate.

fricas
d sin(x)
 (5)
Type: DExpression?
fricas
d acos(x)
 (6)
Type: DExpression?
fricas
d(x^(-1/3))
 (7)
Type: DExpression?

Clearly, our code was designed with the chain rule, so we can also differentiate more complicated composed functions.

fricas
d sin acot exp(x^(1/3))
 (8)
Type: DExpression?

## Comments on the above program

• Comments in SPAD are introduced with double minus -- and reach till the end of the line.
• SPAD knows docstrings. They are started with double plus ++ and reach until the end of line. Docstrings must be attached to actual program code. As can be seen from above, one can document whole categories and domains, as well as function signatures. These docstrings are shown in the HyperDoc system.
• Although MAbelianAdditiveGroup is a category, it has an implementation part (introduced with the keyword add). Function implementations that go in such a place are generic in the sense that they cannot refer to any specific representation of a corresponding domain.
• There is a convention is SPAD that all functions that decide something (like variable?), i.e., return a Boolean value, should end with a question mark. The question mark is part of the name.
• The form
      B => T
E1
E2


      if B then
T
else
E1
E2

• A form like
      for z in l for j in 1..len repeat foo(z, j)


runs in parallel through each element z of the list l and from j=1 until j=len. It will call

      f(l.1, 1)
f(l.2, 3)
f(l.3, 3)
...


The function foo(z,j) is called at most len times, but if the list has fewer than len elements, then there will be only as many invocations as the list has elements.

• The function # is defined for lists in the FriCAS library. It returns the length of a list.
• The error function prints the string that it gets as an argument to the console and then aborts the program. There is not yet an exception mechanism is SPAD.
• In SPAD one can write f n or f.n for f(n). The difference is in how multiple compositions associate. We have:
• x.y.z means (x.y)(z) which is the same as (x(y))(z), and
• x y z means x(y(z)).