From: katherle@daisy.uwaterloo.ca (Kate Atherley)
Subject: Re: Eliminating zeroes from multi-valued expressions in MAPLE
Organization: University of Waterloo
Date: Tue, 17 Nov 1992 17:35:08 GMT
===============================================================================

In article <924816194818@ibm3090.bham.ac.uk> PEARCES@UK.AC.BHAM.IBM3090 writes:
>I`m a MAPLE user with a big problem.
>The following expression has
>zeroes at phi_d := {.9807315053, 2.160861149}.
>
>            162265091            2    446229  1/2
>   (223 (- ----------- sin(phi_d)  + --------)
>               4000                     11
>
>                162265091            2    7689337  1/2
>      + 100 (- ----------- sin(phi_d)  + ---------)   )
>                   4000                      85
>
>      ^ 2
>Does anyone know how to find the function left after
>eliminating each pole. My MAPLE procedure fails to work properly
>even if I only wish to evaluate the value of the expression at each
>pole when the pole has been removed. Simple dividing the expression
>by an (phi_d-0.98073...) and taking the limit of the resultant expression
>as phi_d approaches 0.98703.... Doesn`t work in MAPLE either.

Although Maple's solve() function does return the roots you give (in
exact rational form until you apply evalf()), these are, in fact, not
roots of your function.  They are spurious roots introduced by squaring.

The problem is that the expression goes complex near Pi/2.  This is why
the spurious roots appear.  And where the expression is real, it is on the
order of 10^10.  Thus there are no real roots to your expression, which
explains why such things as limit have difficulty when you try to divide
out the alleged roots near .98 and 2.16.

Here is a maple session showing some of the details:


    |\^/|      MAPLE V
 .._|\|   |/|_.  Copyright (c) 1981-1990 by the University of Waterloo.
 \  MAPLE  /   All rights reserved.  MAPLE is a registered trademark of
 <____ ____>   Waterloo Maple Software.
      |        Type ? for help.
# Increase Digits to protect against cancellation problems
> Digits:=40:

# Define the subexpressions
> p1 := phi_d -> -162265091/4000*sin(phi_d)^2+446229/11;
                                                      2
            p1 := phi_d -> - 162265091/4000 sin(phi_d)  + 446229/11

> p2 := phi_d -> -162265091/4000*sin(phi_d)^2+7689337/85;
                                                      2
            p2 := phi_d -> - 162265091/4000 sin(phi_d)  + 7689337/85

# Here the original problem:
> p := (223*p1^(1/2) + 100*p2^(1/2))^2;
                                     1/2         1/2 2
                         p := (223 p1    + 100 p2   )

# Maple's attempt to find roots of p (but they are not roots)
> solve(p(phi_d));
                            %1, %1, Pi - %1, Pi - %1

                                                              1/2
%1 :=     arccos(1/70912927803729 1556855523643604433146349441   )

> evalf(["]);
                [.9807315051080000303322044017505411924602,

                    .9807315051080000303322044017505411924602,

                    2.160861148481793208130438981528961691737,

                    2.160861148481793208130438981528961691737]

# Get the roots of the first piece
> solve(p1(phi_d));
               180                  1/2              180                  1/2
  - arcsin(---------- 98331022495090   ), arcsin(---------- 98331022495090   )
           1784916001                            1784916001

> evalf(["]);
                 [-1.570772657184545921339569766375277706391,

                     1.570772657184545921339569766375277706391]

# Periodicity gives another (family of) roots
> evalf(Pi-"[2]);
                   1.570819996405247317123073616904225177806

# How far apart are they?
> "-""[2];
                    .000047339220701395783503850528947471415

# Get the roots of the second piece
> solve(p2(phi_d));
                               20                      1/2
                  - arcsin(---------- 42422172913178678   ),
                           2758506547

                                 20                      1/2
                      arcsin(---------- 42422172913178678   )
                             2758506547

# No real roots:
> evalc(evalf(["]));
            [ - 1.570796326794896619231321691639751442099

                  - .9564233764000061305075721869855452404189 I,

                1.570796326794896619231321691639751442099

                     + .9564233764000061305075721869855452404189 I]

> p1(Pi/2);
                                    -1/44000

> evalf(");
                 -.00002272727272727272727272727272727272727273

> p2(Pi/2);
                                   3392963053
                                   ----------
                                      68000

> evalf(");
                   49896.51548529411764705882352941176470588

> p1(Pi);
                                     446229
                                     ------
                                       11

> evalf(");
                   40566.27272727272727272727272727272727273

> p2(Pi);
                                    7689337
                                    -------
                                       85

> evalf(");
                   90462.78823529411764705882352941176470588

> p(Pi/2);
               / 223         1/2                   1/2      1/2\2
               |----- I 44000    + 1/680 3392963053    68000   |
               \44000                                          /

# This is about as small as the function gets
> evalc(evalf("));
                                                           9
               .4989651537227366310160427807486631016040*10

                    + 47494.53383446247774655548050782504198492 I

> p(Pi);
                 /223       1/2   1/2    20         1/2   1/2\2
                 |--- 446229    11    + ---- 7689337    85   |
                 \ 11                    17                  /

> evalf(");
                                                             10
                 .5623743076386678984956221653763304287514*10

> quit

Kate (with thanks to Dave Hare and Greg Fee)
Waterloo Maple Software Technical Support
support@maplesoft.on.ca
First of all may I express my utmost thanks to Kate et al.
Unfortunately the above does not help me with my problem that much.
Therefore I will try to explain it more fully.
Sadly this will mean using terms and concepts with which I am still
as yet not fully comfortable.

Firstly my original expression contains square roots.
These are multi-valued and hence we have 2 Rieman sheets for the
rooted parts of the expression.
Since there are 2 of these does this mean that the expression without
the ^2 has 4 Rieman sheets and after squaring this expression
has 2?
I have worked out by hand that the `spurious' roots quoted above
correspond to the expression when we choose the negative root
for one of the square roots and the positive for the other.
This would mean that the `spurious' roots are on a different Rieman
sheet to that normally evaluated by MAPLE.
It is this situation I am most interested in since it has physical
significance in the problem I am tackling.

So the question in a nut shell is how to change the original expression
so that it corresponds to the `spurious' roots and then eliminate the
zero to find the remaining function (or just its value at the zero
if this is not possible).
Since this is to form part of a problem solving system
I need a routine to do this automatically.
Any further help would be greatly appreciated.
Yours Gratefully,
Stephen Pearce.
P.S. I am using MAPLE 4.2 on an IBM3090 and so did not fully understand
the MAPLE V syntax in the last message.
Newsgroups: sci.math.symbolic