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
Newsgroups: sci.math.symbolic