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

[Sumant S.R. Oemrawsingh]

Hi guys,

I recently had some time to try some axiom stuff again. I've tried to get
the piecewise functions to work.

First thing I did, was try with a heaviside function:

(1) -> h(k,x)==1/2 + 1/%pi * atan(k*x)
h(k,x)==1/2 + 1/%pi * atan(k*x)
                                                                   Type: Void

Now, in the limit of k->infinity, this is exactly the heaviside function. As
such, it is possible to write a piecewise function using this:

(2) -> f(k,x)==x**2 * (h(k,x)-h(k,-x))
f(k,x)==x**2 * (h(k,x)-h(k,-x))
                                                                   Type: Void

which gives exactly the results that are expected:

(3) -> draw(limit(f(k,x),k=%plusInfinity),x=-5..5)
draw(limit(f(k,x),k=%plusInfinity),x=-5..5)
   Compiling function h with type (Variable k,Variable x) -> Expression
      Integer 
   Compiling function h with type (Variable k,Polynomial Integer) -> 
      Expression Integer 
   Compiling function f with type (Variable k,Variable x) -> Expression
      Integer 
   Compiling function %A with type DoubleFloat -> DoubleFloat 
   Graph data being transmitted to the viewport manager...
   AXIOM2D data being transmitted to the viewport manager...

   (3)  TwoDimensionalViewport: "x*DABS(x)"
                                                 Type: TwoDimensionalViewport


However, integration still doesn't work, which is rather odd:

(4) -> g(x)==limit(f(k,x),k=%plusInfinity)
g(x)==limit(f(k,x),k=%plusInfinity)
                                                                   Type: Void
(5) -> integrate(g(x),x=-1..1)
integrate(g(x),x=-1..1)
   Compiling function g with type Variable x -> Union(OrderedCompletion
      Expression Integer,Record(leftHandLimit: Union(OrderedCompletion 
      Expression Integer,"failed"),rightHandLimit: Union(
      OrderedCompletion Expression Integer,"failed")),"failed") 

   (5)  potentialPole
                                         Type: Union(pole: potentialPole,...)
(6) -> integrate(g(x),x=-1..1,"noPole")
integrate(g(x),x=-1..1,"noPole")
(6) -> 
   (6)  "failed"
                                                Type: Union(fail: failed,...)

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. Naturally, the
integral over the heaviside function itself also doesn't work, as there is a
pole in the path of integration at x=0.

Now, since the heaviside function is special in a way, that we can remove
the function from the integration, as long as we adjust the boundaries
appropriately, I would like to somehow implement this. Or at least,
integration over a piecewise defined function should automatically perform
this adjustment of the integration boundaries.

Closely linked to that is integration over the function abs(x). Apparently,
there is a potential pole there as well. In a sense, abs(x) can be written
using heaviside functions as well (or the other way around, as I did
previously). Thus, integration over a function like abs(x) should work
properly by splitting the integration into two separate integrals, when x=0
is in the integration domain.

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.

I also have difficulty finding some sort of example of how to do these types
of things in spad; the book and other materials I found on the web don't
really explain how to extend the algebraic functionality, or I haven't
looked in the proper places... So any pointers (web-wise or source-wise)
would be greatly appreciated. If I can calculate the integral of such
simple, piecewise functions by hand, I find that axiom should be able to do
so as well (symbolically).

Regards,

Sumant




On Tue, Jun 05, 2007 at 08:49:56PM +0200, Sumant S.R. Oemrawsingh wrote:
> Hi Bill,
> 
> Thanks for your elaborate and clarifying answers. I agree with you that I
> will probably have to help build what I want. I don't have much time to
> spend on it, but I am looking at it currently and will hopefully be
> contributing something in time.
> 
> Regards,
> 
> Sumant
> 
> On Sun, Jun 03, 2007 at 02:03:50AM -0400, Bill Page wrote:
> > On June 2, 2007 10:21 PM Sumant S.R. Oemrawsingh wrote:
> > > 
> > > Thanks for your answer. It is now more clear to me what 
> > > happens. As I am completely new to Axiom, I probably don't
> > > understand the full philosophy behind it yet, so I wish to
> > > ask a bit more on this, if you don't mind.
> > 
> > That is the purpose of this email list. Everyone is free to
> > ask and to answer questions here but if we discuss philosophy
> > I am afraid not everyone will agree. I believe that Axiom was
> > designed with a certain overall philosophy but that it was
> > not necessarily always applied consistently or in the same
> > manner by different people. For sure not everyone is going
> > to agree with my personal opinions about exact what this
> > "Axiom philosophy" is. I think that to properly discuss some
> > issues of philosophy requires a deep and thorough understanding
> > of the subject but not everyone has the patience or skill to
> > work at this level. Many people have very pragmatic reasons
> > for wanting to use a computer for doing mathematics - they have
> > some job they want to do and do not care too much about the
> > subject from an abstract "computer science" point of view. But
> > in my opinion the designers of Axiom did try very hard to take
> > the more general view while the designers of most other
> > superficially similar systems (e.g. reduce, maxima, maple,
> > mathematica, etc.) did not.
> > 
> > > 
> > > > You said you expected the 'f(x)' would represent some kind
> > > > of expression defined in a piecewise manner, but Axiom does
> > > > not currently have any expressions of this kind.
> > > 
> > > This is very strange, since in the Axiom book, there is a 
> > > whole section on piece-wise functions. Why have them, if
> > > they are not supported, so to say?
> > 
> > You have a very good point. Yes Axiom does have piece-wise
> > defined *functions* but it does not have piece-wise *expressions*.
> > It is important to make the distinction between function and
> > expression. In fact going deeper one can see that Axiom uses
> > the concept of "expression" in two different ways. In my
> > opinion understanding this is also an important part of the
> > Axiom philosophy.
> > 
> > For example in Axiom you can write the expression:
> > 
> >   (1) -> if x>0 then 1 else 0
> > 
> >     (1)  1
> >                               Type: PositiveInteger
> > 
> > The result of evaluating this expression depends on the value
> > assigned to x. This is a "piece-wise" expression in the sense
> > we have been discussing *BUT* this is a part of the Axiom
> > programming language - it is not part of any mathematical
> > object currently implemented in the Axiom library. The Axiom
> > library includes things like polynomials and even the domain
> > 'Expression' containing expressions like 'x+1' or 'sin(1)'
> > but there is currently no domain in which we can find an
> > expression like 'if x>0 then 1 else 0'.
> > 
> > > If it can be defined and yields correct answers when called 
> > > with integers or even expressions yielding specific values
> > > (say, sqrt(5) or cos(1/2)), why can't it be integrated over
> > > a domain where the variable takes on specific values for
> > > which all the definitions are required to be evaluated?
> > > 
> > 
> > Do you understand the difference between writing?
> > 
> >   draw(sin(x),x=-5..5)
> > 
> > and
> > 
> >   draw(sin, -5..5)
> > 
> > Martin showed how you could use draw in this way to display a
> > piecewise defined function. Unfortunately Axiom's symbol
> > integration operator does not currently allow a function as
> > argument but Martin did gave an example of using 'Romberg' to
> > do numerical integration of functions.
> > 
> > In principle it would be possible to extend Axiom's integration
> > and differentiation operators to work symbolically on functions
> > but this is a non-trivial problem. I believe one Axiom developer
> > (Gabriel Dos Reis one of his students) has done some work on
> > symbol differentiation of functions.
> > 
> > It is also possible in principle to deal with piece-wise
> > expressions as mathematical objects in Axiom. Tim Daly has
> > discussed the development of a general "proviso" feature in
> > Axiom that would make this possible. But even fairly extensions
> > of the existing 'Expression' domain could probably handle
> > simple cases of piece-wise expressions of the sort you have
> > discussed so far.
> > 
> > > Thanks to your explanation, I understand why x < 0 (Variable, 
> > > Polynomial Integer, Symbol,...) *does not* return true, but
> > > I do not understand why it *should* not return true (from a
> > > design point of view).
> > 
> > Perhaps you are objecting to the decision of the Axiom library
> > developers to define '<' in terms of lexical ordering for
> > polynomials. Although this is natural in a certain context,
> > I agree that that was probably a poor choice on their part.
> > 
> > > But from a mathematical point of view (imho), the expression
> > > x < 0 should not be false; it should be true by default,
> > > unless it has been specified that x is an element of some
> > > set, since if it is *not* defined, it can have *any* value,
> > > also negative. Therefore, when one expects to test the
> > > domain of x, such as in an expression like this, it should
> > > be true by default.
> > 
> > I strongly disagree with this. Most people would read 'x<0'
> > as a logical predicate. It is true or false depending on
> > the value and/or quantification associated with x. You are
> > right to be concerned about the domain of x but unfortunately
> > Axiom currently has no way to evaluate logical predicates of
> > that kind.
> > 
> > > 
> > > So, I wish to tell axiom that x is an element of, say, R.
> > > I understand that I should declare the type, in my example
> > > Float -> Float, like so:
> > > 
> > > (1) -> f: (Float) -> Float
> > >                                                               
> > >      Type: Void
> > > (2) -> f(x|x<=0)==-x**2   
> > >                                                               
> > >      Type: Void
> > > (3) -> f(x|x>0)==x**2     
> > >                                                               
> > >      Type: Void
> > > (4) -> f(5)
> > >    Compiling function f with type Float -> Float 
> > > 
> > >    (4)  - 25.0
> > > 
> > > (5) -> f(5**(1/3))
> > 
> > >    (5)  2.9240177382 128660655
> > >                                                               
> > >     Type: Float
> > >
> > 
> > Right.
> >  
> > > How do I get a symbolic answer now? I don't care what the 
> > > answer up to n digits is, I want to know the actual, exact
> > > answer.
> > 
> > But you specified the domain 'Float' that is a domain where
> > numbers are represented by an approximation. If you want
> > symbolic answers you need to specify a different domain,
> > e.g. Expression Integer.
> > 
> > > This is even worse:
> > > 
> > > (6) -> draw(f(x),x=-1..1)
> > >    Conversion failed in the compiled user function f .
> > >  
> > >    Cannot convert from type Symbol to Float for value
> > >    x
> > > 
> > 
> > Write it this way:
> > 
> >   draw(f,-1..1)
> > 
> > When you write 'f(x)' you are asking for a piece-wise
> > expression but that is not currently possible in Axiom. But
> > if you operate directly on the piece-wise function then you
> > can get what you want.
> > 
> > > (6) -> integrate(f(x),x=-1..1)
> > >    Conversion failed in the compiled user function f .
> > >  
> > >    Cannot convert from type Symbol to Float for value
> > >    x
> > >
> > 
> > Unfortunately you can not currently write:
> > 
> >   integrate(f, -1..1)
> > 
> > except in the case where you would like to do numerical
> > integration (e.g. Romberg).
> >  
> > > Granted, I really don't care that much about drawing the 
> > > function. It is more of a visual aid to see that I defined
> > > my function correctly. For the example, this is not really
> > > necessary, but some of the functions I wish to perform
> > > calculations with, have some wild definitions, with many
> > > pieces.
> > 
> > For drawing that is ok.
> > 
> > > However, integration is very important, since that is what
> > > I really want to do in the end. Even if it manages to give
> > > me an answer of the integral, I still don't want to know
> > > the answer estimated to n digits, I would like the exact
> > > answer.
> > >
> > 
> > I think that exact integration of piece-wise functions in
> > general is a rather difficult problem.
> >  
> > > So the real question I would like to ask is, how can I
> > > tell axiom that, during function definition, my variable
> > > belongs to R (or C), while all operations still remain
> > > symbolic?
> > 
> > Right now there is no way to do this in Axiom but that is
> > (more or less) what Tim Daly means by a 'proviso'. This is
> > still a research topic in Axiom.
> > 
> > > How can I get it to first evaluate the value of x, figure
> > > out which function to use, and then return the value
> > > symbolically again?
> > > 
> > 
> > I don't know.
> > 
> > > I have tried by using "numeric(x) > 0" in the conditions,
> > > but then integration won't work since numeric won't work on
> > > symbols etc. Thanks to your explanation, I already expected
> > > that before I even tried it. On a side note, f(%pi) will
> > > work correctly then.
> > > 
> > > Which leads me to the following: why is %e fundamentally 
> > > different from %pi? Or, since Martin sort of answered that
> > > question already, I should ask: why *should* %e be
> > > fundamentally different from %pi? 
> > 
> > I don't think it should be. This is just an accident of a
> > series of separate design decisions. As Martin said, this
> > should probably be corrected by a more or less complete
> > redesign of how Axiom handles constants. This would amount
> > to another Axiom research project.
> > 
> > > 
> > > But I digress. The most important question is, how to define
> > > a function such that it first evaluates the value of x, then
> > > chooses which case to use and then perform the action?
> > > 
> > > In integrate(), I would expect it would look at the limits,
> > > perform an indefinite integral of the relevant cases, and
> > > then fill in the limits for each case correctly.
> > 
> > Perhaps you could describe this algorithm in more detail.
> > Can you write some Axiom scripts that perform this sort
> > of computation?
> > 
> > > For draw, I would expect something similar.
> > > 
> > > I've been doing some other stuff with axiom (no problems 
> > > there!), and think it is really amazing. That is also the
> > > reason that I find this particular point so frustrating. So
> > > if you or any one could help, I would be very much obliged.
> > > 
> > 
> > I hope this reply is of some help even though for the most
> > part it is saying that what you want from Axiom you will
> > probably have to help build for yourself. :-)
> > 
> > Regards,
> > Bill Page
> > 
> > 
> 
> -- 
> Sumant S. R. Oemrawsingh
> address@hidden
> 
> 
> _______________________________________________
> Axiom-mail mailing list
> address@hidden
> http://lists.nongnu.org/mailman/listinfo/axiom-mail

-- 
Sumant S. R. Oemrawsingh
address@hidden