Symbolic Mathematical Computing and
Problem Solving Environments
Richard Fateman
Computer Science Division, EECS
University of California, Berkeley
February, 1994
Outline
1. Introduction
Symbolic > Numeric
2. Symbolic Mathematics Components
Can we extract them?
Program that manipulate programs
Derivatives, integrals
Closed form solutions
Semi-symbolic solutions
Exact, high-precision, interval ... arith.
Code generation for special arch.
Support for proofs, etc
Documentation/Output
Databases
Specific toolboxes: FEA, etc.
3. Symbolic Manipulation Systems as Glue
An argument for Lisp as scripting
and symbolic language
4. Objectives for Symbolic Computing
5. Future tools and problem areas
INTRODUCTION
What is Symbolic Computation and how is it different?
Symbolic includes Numeric, Graphical, ...
Computer Algebra systems are becoming
interactive environments
(Mathematica, Maple, Macsyma, Axiom)
Computer Algebra systems are being joined with
other environments (e.g. MathCAD+Maple)
PART I: SYMBOLIC COMPONENTS
Toolkit approach is difficult:
It is hard to tease out individual pieces
[storage model, abstraction, run-time semantics]
(and futile to do so, in some cases)
{Commercial vs. Research objectives}
-------------------
Program that manipulate programs
* Assemblers, interpreters,
(pre-)Compilers, macro-expansion, etc.
The symbolic view is that:
expressions <==> programs <==> data
Consider this section
of a Bessel Function evaluation pgm. (from Numerical Recipes)
...
DATA Q1,Q2,Q3,Q4,Q5,Q6,Q7,Q8,Q9/0.39894228D0,-0.3988024D-1,
* -0.362018D-2,0.163801D-2,-0.1031555D-1,0.2282967D-1,
* -0.2895312D-1,0.1787654D-1,-0.420059D-2/
...
BESSI1=(EXP(AX)/SQRT(AX))*(Q1+Y*(Q2+Y*(Q3+Y*(Q4+
* Y*(Q5+Y*(Q6+Y*(Q7+Y*(Q8+Y*Q9))))))))
...
Bessel Function/ Polynomial Eval in Lisp
(program --> better program )
(setf
bessi0
(* (/ (exp ax) (sqrt ax))
(poly-eval y (0.39894228d0 -0.3988024d-1 -0.362018d-2 0.163801d-2
-0.1031555d-1 0.2282967d-1 -0.2895312d-1 0.1787654d-1
-0.420059d-2))))
extract the polynomial part, and reformulate
(let* ((z (+ (* (+ x -0.447420246891662d0) x) 0.5555574445841143d0))
(w (+ (* (+ x -2.180440363165497d0) z) 1.759291809106734d0)))
(* (+ (* x (+ (* x (+ (* w (+ -1.745986568814345d0 w z))
1.213871280862968d0))
9.4939615625424d0))
-94.9729157094598d0)
-0.00420059d0))
uses 6 multiplies, 9 adds, not 8 of each (Better than just
Horner's rule!)
Euler equation
(expression --> program(s))
The Euler equation is a favorite benchmark of Celestial Mechanics
symbolic calculation programs.
Systems without Poisson series do poorly on it.
E = u + e sin (E)
as commonly solved iteratively (for small e) gives a
4th order expansion for E= u+A as:
4 3 2 4
e sin(4 U) 3 e sin(3 U) (12 e - 4 e ) sin(2 U)
A = ----------- + ------------- + -----------------------
4 3 8 24
3
(24 e - 3 e ) sin(U)
+ --------------------
24
Euler equation (continued 1)
In Fortran as rendered by Mathematica:
FortranForm=
- (24*e - 3*e**3)*Sin(U)/24 + (12*e**2 - 4*e**4)*Sin(2*U)/24 +
- 3*e**3*Sin(3*U)/8 + e**4*Sin(4*U)/3
(In general, Mathematica is outright wrong in conversion
into Fortran since it produces forms like 1/3 )
Call... Expand[N[%]]
FortranForm=
- 1.*e*Sin(U) - 0.125*e**3*Sin(U) + 0.5*e**2*Sin(2.*U) -
- 0.1666666666666666*e**4*Sin(2.*U) + 0.375*e**3*Sin(3.*U) +
- 0.3333333333333333*e**4*Sin(4.*U)
Apparently Mathematica's Fortran has some other errors: it
produces text with multiplications by 1., and has decided exactly
how many 3's are relevant for your Fortran compiler.
Euler equation (continued 2)
Maple produces
t0 = E**4*sin(4*U)/3+3.0/8.0*E**3*sin(3*U)+(12*E**2-4*E**4)*sin(2*
#U)/24+(24*E-3*E**3)*sin(U)/24
or after floating-point conversion using evalf
t0 = 0.3333333E0*e**4*sin(4.0*U)+0.375E0*e**3*sin(3.0*U)+0.4166667
#E-1*(12.0*e**2-4.0*e**4)*sin(2.0*U)+0.4166667E-1*(24.0*e-3.0*e**3)
#*sin(U)
Maple has also decide how many 3's in 1/3, though this can be changed.
After rearranging using Horner's rule (convert(expression,horner,[e])}))
t0 = (sin(U)+(sin(2*U)/2+(3.0/8.0*sin(3*U)-sin(U)/8+(-sin(2*U)/6+s
#in(4*U)/3)*e)*e)*e)*e
Note sin(u) occurs twice..
Euler equation (continued 3)
Macsyma can do about the same but also can pick out common
subexpressions. But notice that
s1 := sin(u)
c1 := cos(u)
s2 := 2*s*c
c2 := 2*c*c-1
s3 := s*(2*c2+1)
s4 := 2*s2*c2
gives s2 = sin(2 u), s3 = sin(3 u), s4 = sin(4 u) ?
None of the computer systems notice this.
Derivatives, Integrals
Careful approaches (Griewank and Corliss) to derivatives can help a
great deal more than Computer Algebra systems. Naively sticking in
symbolic derivatives is often a bad idea.
Simple example: Evaluating a polynomial AND its derivative can be done
a lot faster by a modest change to Horner's reccurence than by using a
new expression.
Integrals
Closed form or quadrature? which is faster?
integrate(1/(1+z^5),z), or even worse, integrate(1/(1+z^64),z)
can be integrated, but the form is horrendous.
Whose closed form do we use
Integrate[1/x,{x,a,b}] or
integrate(1/x,x,a,b); or
int(1/x,x=a..b);
all give log(b)-log(a).
Not log(b/a) or if 1/2** a sure-fire tool.
Closed form solutions
The closed form symbolic expression, wearing singularities
on its face: The exactly zero expression. Proving identities.
Semi-symbolic solutions
Taylor, Laurent, Fourier, Asymptotic- series, perturbation
expansion; expansions in terms of sets of orthogonal polynomials
Exact, high-precision, interval ... arith.
Arithmetic on objects of variable size is offered:
Most systems support exact integer and rational computing.
This breaks down for exponential, log, trignometric function
computing --> arbitrary-precision floats.
exp(pi*sqrt(163))=262537412640768743.9999999999992500726.
New algorithms can be developed, e.g. for root-finding.
Old analyses for fixed precision are no longer valid.
Code generation for special architectures
Support for derivations, proofs, reasoning,
Justifications, etc.
Documentation:
Output as TeX
Notebooks (Mathematica, Maple)
Spreadsheets (Theorist, MathCAD)
Graphics into AVS, other graphics packages
Specific Toolboxes e.g.
Finite Element Analysis
PART III
Symbolic Manipulation Systems as Glue
Fortran or C cannot "naturally" talk about objects
the way (for example) Lisp or Axiom can.
Matrix languages (Matlab etc) are also deficient.
Character strings in pipes or RPC don't work
well enough.
Is there a role for Lisp?
Lisp as scripting and symbolic language
Other scripting languages.
PART IV
Future objectives for Symbolic Tools and Environments
* Integration of data bases of properties, formulas, algorithms.
* Refined application packages.
* Geometric, constraint-based reasoning.
* Linear and non-linear inequalities.
* Very fast constructive mathematics algorithms using all the tools
at our disposal:
floating point, parallel/distributed techniques (e.g.
Strand-88 and Linda- Maple)
* Management of queries into large data bases (e.g. algebraic
or logical simplification of queries)
* Access to networked data and algorithms (e.g. the world's
best integration program; ODE program generator; closed-form
summation etc.)
*****IMPORTANT TEST CASES*******
Future Objectives for Glue
* Ease of use -- e.g. AVS or Khoros style box programming.
* Ease of use -- clean functional /scription language.
* Enhances portability
* Access to "everything" in other languages.
Other Symbolic Math Activities
* OCR, Parsing, storage and retrieval of mathematical
formulas (largish database for matching)
* Theorem proving (plane geometry)
* Common Lisp as source language for numerical
computing, retrospective diagnostics, etc.
* High-speed Parallel polynomial arithmetic.
Novel data structures for serial or parallel
manipulation.
* Honest plotting (guarantees via interval arithmetic
that what you see is not misleading)
**