Re: [fricas-devel] Re: [Axiom-math] Selection of roots

[Bill Page]

On Mon, Nov 3, 2008 at 7:55 PM, Alejandro Jakubi wrote:
> Billl Page wrote:
>> One thing that would be very useful is if someone like you were able
>> to contribute some documentation like this that could be added to
>> the axiom-wiki and perhaps also the Axiom book.
>
> May be. But wait. I am slow and it takes me a very long time to
> understand these kind of things. At present I do not understand even
> the basics of the Axiom system. And I do not find sensible that try
> contributing to this documentation until I beleive that I understand the
> subject enough.
>

There are different kinds of documentation. I expect that the best
time to write a tutorial (for other people at a similar stage) is
while one is still learning, or at least while it is still "fresh".

> You see, it is over 15 years that I am using Maple, and it is only now
> that I feel that I understand enough about a few subjects to begin
> writing on them in the Maple wiki
> (http://mapleadvisor.com/cgi-bin/moin.cgi/StartingPage)
>

Well, Tim Daly states that the original Axiom project has a 30 year
time horizon, so I guess you still have time. :-)

>> BTW, what is you basic source of information about Axiom?
>
> I think that there is nothing that I could call a basic source of
> information about Axiom.
>
> For several years I have been collecting all sort of documentation,
> including the Axiom book (in different versions, formats, etc),
> tutorials, posts to the lists, wiki pages, etc.
>
> And now I have began trying to "decode" the output of commands
> like ')show' (without much success yet).
>

Hyperdoc is also a very important source of information. It is often
easier and faster to find what you are looking for there. But
ultimately (for better or worse, mostly worse) the last resort is the
algebra source code itself.

> In most cases, none of them by itself, provides an answer to my
> howto questions. Hence, currently, I am collecting such dispersed
> bits of information plus examples in TeXmacs sessions to help orient
> myself in this search.
>

Great. I think you work would also be of some interst to others.

>> I believe that all you are really interested in are the roots of the
>> numerator of the expression (provided of course that they are not
>> also roots of the denominator), right?
>
> Right.
>
>> (3) -> S:=allRootsOf(numer fp)$RealClosure(Fraction Integer)
>>
>>  (3)  [%A9,%A10]
>>                                     Type: List RealClosure Fraction
>> Integer
>> (4) -> approximate(S.1,1/10^20)::Float
>>
>>  (4)  - 6.7957899636 620037966
>>                                                                 Type:
>> Float
>
> The final step in these examples would be to recover the algebraic
> expression of these selected roots (instead of a float evaluation).
>

The documentation says RealClosure characterizes solutions uniquely
bound by intervals open on the right, i.e. just one real-valued
solution is known to lie at or greater than the "left" value and less
than the "right" value. It does not actually find exact solutions.
'radicalSolve' on the other hand finds exact symbolic values (up to
4th order) but these can be difficult to evaluate.

It is sometimes easy to find the exact real solutions directly:

  A:=select(x+->retractIfCan(complexNumeric rhs x) case Float,radicalSolve(fp));

such as in this particular example, but this method could silently
skip some roots due to rounding errors.

Here is one possible script that uses RealClosure's interval
characterization to locate the associated symbolic (exact) roots.

---
f:=(x^3+x^2-4*x-4)/(2*x^2+7*x-4)
fp:=differentiate(f,x)
-- locate real roots
S:=allRootsOf(numer fp)$RealClosure(Fraction Integer);
-- find exact solutions
R:=radicalSolve(fp);
-- test that s is bound by s
bound?(x,s) == (a:=complexNumeric rhs x; _
                   imag a < 10^-digits() and _
                   real a >= left(mainCharacterization s)::Float and _
                   real a < right(mainCharacterization s)::Float)
-- find exact real roots
A:=[  (B:=select(x+->bound?(x,s),R); _
          #B=1 => B.1; _
          error "failed unexpectedly") _
        for s in S ];
-- evaluate them
map(x+->real complexNumeric rhs x, A)
---

I haven't tested this extensively. I suppose there are might be some
cases where it could also fail due to floating point rounding errors,
but it seems to work for the two examples you have given so far.

Regards,
Bill Page.