Re: [Axiom-developer] Re: [Axiom-mail] Defining piece-wise functions and drawing, integrating, ...

[Waldek Hebisch]

Martin Rubey wrote:
> Dear Sumant,
> 
> "Sumant S.R. Oemrawsingh" <address@hidden> writes:
> 
> 
> > (1) -> h(k,x)==1/2 + 1/%pi * atan(k*x)
> 
> > (2) -> f(k,x)==x**2 * (h(k,x)-h(k,-x))
> 
> > (4) -> g(x)==limit(f(k,x),k=%plusInfinity)
> 
> > What I would think is that, since g(x) is a sum of two functions, the
> > integral would split up into a sum of integrals as well.
> 
> Hm, why should g(x) be a sum of two functions? My axiom gives me upon issueing
> g x
> 
>  x abs(x)
> 
> > (6) -> integrate(g(x),x=-1..1,"noPole")
> > integrate(g(x),x=-1..1,"noPole")
> > (6) -> 
> >    (6)  "failed"
> >                                                 Type: Union(fail: 
> > failed,...)
> 
> This just says that axiom couldn't do the integral.  In particular, bugs 
> aside,
> it is a "proof" that it not expressible as an elementary function.  (The Risch
> algorithm is more or less completely implemented in axiom, and axiom will spit
> out an error message if you hit an unimplemented branch.)
> 

Note, that from Axiom point of view x*abs(x) is _not_ an elementary
function, so Risch algorithm immediatly gives up.

> > So I've been looking a bit into how this could be done in spad. But I've not
> > been able to understand where and how the functionality of such special
> > functions is or can be implemented, or if I somehow would have to extend the
> > definition or workings of the integrate function itself.
> 
> For a start, I'd think it is easier to split up integration boundaries
> semi-automatically, as I indicated in a previous email.
> 
> If you are really interested in the internals, I think that intpm.spad is for
> you, it tells axiom some pattern matching rules.  In my opinion, we should go
> the mupad way here, the source of axiom is a horror, while the one of mupad is
> extremely easy to grasp and to extend.
> 

I do not understand the comment about mupad.  I did not see mupad code,
but source of Axiom integrator is very clear: it implements Risch
algorithm very much as described in research articles.  Basically
the only extra difficulty is due to use of overloading: somethimes
one have to spend some time to discover domain from which given
operation comes.  The algorithm (and expecially math behind it) is
complex, so anybody who does not know the algorithm may easily get
lost.  But given complexity of math I doubt that that it can be
written much clearer.

Concerning integration of piecewise functions: there is an
algorithm which handles a buch of interesting cases. See:

D. J. Jeffrey,  G. Labahn,  M. v. Mohrenschildt,  A. D. Rich,
Integration of the signum, piecewise and related Functions
http://citeseer.ist.psu.edu/jeffrey97integration.html

and

D. J. Jeffrey,  A. D. Rich,   Recursive integration  of
piecewise-continuous functions
http://citeseer.ist.psu.edu/jeffrey97recursive.html

If one wants to implement something similar in Axiom the
first thing to do is to define a domain of piecewise
continuous functions.  Why?  Normal Axiom expressions form
a field, piecewise continuous functions have zero divisors,
so they only form a ring.  Zero divisors, if unhandled
will cause various nonsense results.

In general I belive that people who want to change Axiom
should quickly define their own category/domain/package
and work there.  Once such code is mature we will try
to integrate it.  It would be good to have some
examples for such extensions.  While defining new domains
looks easy one can hit problems very quickly.  For
example, a little package that tries to model adding
special functions:

)abbrev package EPAK MiscFunctions
MiscFunctions() : Exports == Implementation where
  F == Expression Integer
  Exports ==> with
     HypergeomF : (List F, List F, F) -> F
  Implementation ==> add
     ophypergeom := operator("HypergomF"::Symbol)$CommonOperators
     oplist := operator("%list"::Symbol)$CommonOperators
     HypergeomF(a, b, z) == ophypergeom [oplist a, oplist b, z]


One problem is that Axiom operats want expressions as arguments.
But here is natural to allow lists -- othewise one would be
forced to put arbitrary limit on the number of parameters and
define bunch of functions which differ only in the number
of paramters they accept.  Above I "solved" this using
atificial "%list" operator, but this has a drawback that
when I type:

HypergeomF([0, 1, 2], [2, 3], 3)

I get back:

HypergomF(%list(0,1,2),%list(2,3),3)

(and I would prefer list in square brackets).  It also seems
that this definitions will cause problems with derivatives
with respect to parameters...

-- 
                              Waldek Hebisch
address@hidden