home *** CD-ROM | disk | FTP | other *** search
/ NetNews Usenet Archive 1992 #27 / NN_1992_27.iso / spool / sci / math / symbolic / 3003 < prev    next >
Encoding:
Text File  |  1992-11-17  |  6.3 KB  |  179 lines
  1. Newsgroups: sci.math.symbolic
  2. Path: sparky!uunet!charon.amdahl.com!pacbell.com!sgiblab!swrinde!cs.utexas.edu!torn!watserv2.uwaterloo.ca!watdragon.uwaterloo.ca!daisy.uwaterloo.ca!katherle
  3. From: katherle@daisy.uwaterloo.ca (Kate Atherley)
  4. Subject: Re: Eliminating zeroes from multi-valued expressions in MAPLE
  5. Message-ID: <BxvFIK.Lup@watdragon.uwaterloo.ca>
  6. Sender: news@watdragon.uwaterloo.ca (USENET News System)
  7. Organization: University of Waterloo
  8. References: <924816194818@ibm3090.bham.ac.uk>
  9. Date: Tue, 17 Nov 1992 17:35:08 GMT
  10. Lines: 167
  11.  
  12. In article <924816194818@ibm3090.bham.ac.uk> PEARCES@UK.AC.BHAM.IBM3090 writes:
  13. >I`m a MAPLE user with a big problem.
  14. >The following expression has
  15. >zeroes at phi_d := {.9807315053, 2.160861149}.
  16. >
  17. >            162265091            2    446229  1/2
  18. >   (223 (- ----------- sin(phi_d)  + --------)
  19. >               4000                     11
  20. >
  21. >                162265091            2    7689337  1/2
  22. >      + 100 (- ----------- sin(phi_d)  + ---------)   )
  23. >                   4000                      85
  24. >
  25. >      ^ 2
  26. >Does anyone know how to find the function left after
  27. >eliminating each pole. My MAPLE procedure fails to work properly
  28. >even if I only wish to evaluate the value of the expression at each
  29. >pole when the pole has been removed. Simple dividing the expression
  30. >by an (phi_d-0.98073...) and taking the limit of the resultant expression
  31. >as phi_d approaches 0.98703.... Doesn`t work in MAPLE either.
  32.  
  33. Although Maple's solve() function does return the roots you give (in 
  34. exact rational form until you apply evalf()), these are, in fact, not
  35. roots of your function.  They are spurious roots introduced by squaring.
  36.  
  37. The problem is that the expression goes complex near Pi/2.  This is why
  38. the spurious roots appear.  And where the expression is real, it is on the
  39. order of 10^10.  Thus there are no real roots to your expression, which
  40. explains why such things as limit have difficulty when you try to divide
  41. out the alleged roots near .98 and 2.16.
  42.  
  43. Here is a maple session showing some of the details:
  44.  
  45.  
  46.     |\^/|      MAPLE V
  47. ._|\|   |/|_.  Copyright (c) 1981-1990 by the University of Waterloo.
  48.  \  MAPLE  /   All rights reserved.  MAPLE is a registered trademark of
  49.  <____ ____>   Waterloo Maple Software.
  50.       |        Type ? for help.
  51. # Increase Digits to protect against cancellation problems
  52. > Digits:=40:
  53.  
  54. # Define the subexpressions
  55. > p1 := phi_d -> -162265091/4000*sin(phi_d)^2+446229/11;
  56.                                                       2
  57.             p1 := phi_d -> - 162265091/4000 sin(phi_d)  + 446229/11
  58.  
  59. > p2 := phi_d -> -162265091/4000*sin(phi_d)^2+7689337/85;
  60.                                                       2
  61.             p2 := phi_d -> - 162265091/4000 sin(phi_d)  + 7689337/85
  62.  
  63. # Here the original problem:
  64. > p := (223*p1^(1/2) + 100*p2^(1/2))^2;
  65.                                      1/2         1/2 2
  66.                          p := (223 p1    + 100 p2   )
  67.  
  68. # Maple's attempt to find roots of p (but they are not roots)
  69. > solve(p(phi_d));
  70.                             %1, %1, Pi - %1, Pi - %1
  71.  
  72.                                                               1/2
  73. %1 :=     arccos(1/70912927803729 1556855523643604433146349441   )
  74.  
  75. > evalf(["]);
  76.                 [.9807315051080000303322044017505411924602,
  77.  
  78.                     .9807315051080000303322044017505411924602,
  79.  
  80.                     2.160861148481793208130438981528961691737,
  81.  
  82.                     2.160861148481793208130438981528961691737]
  83.  
  84. # Get the roots of the first piece
  85. > solve(p1(phi_d));
  86.                180                  1/2              180                  1/2
  87.   - arcsin(---------- 98331022495090   ), arcsin(---------- 98331022495090   )
  88.            1784916001                            1784916001
  89.  
  90. > evalf(["]);
  91.                  [-1.570772657184545921339569766375277706391,
  92.  
  93.                      1.570772657184545921339569766375277706391]
  94.  
  95. # Periodicity gives another (family of) roots
  96. > evalf(Pi-"[2]);
  97.                    1.570819996405247317123073616904225177806
  98.  
  99. # How far apart are they?
  100. > "-""[2];
  101.                     .000047339220701395783503850528947471415
  102.  
  103. # Get the roots of the second piece
  104. > solve(p2(phi_d));
  105.                                20                      1/2
  106.                   - arcsin(---------- 42422172913178678   ),
  107.                            2758506547
  108.  
  109.                                  20                      1/2
  110.                       arcsin(---------- 42422172913178678   )
  111.                              2758506547
  112.  
  113. # No real roots:
  114. > evalc(evalf(["]));
  115.             [ - 1.570796326794896619231321691639751442099
  116.  
  117.                   - .9564233764000061305075721869855452404189 I,
  118.  
  119.                 1.570796326794896619231321691639751442099
  120.  
  121.                      + .9564233764000061305075721869855452404189 I]
  122.  
  123. > p1(Pi/2);
  124.                                     -1/44000
  125.  
  126. > evalf(");
  127.                  -.00002272727272727272727272727272727272727273
  128.  
  129. > p2(Pi/2);
  130.                                    3392963053
  131.                                    ----------
  132.                                       68000
  133.  
  134. > evalf(");
  135.                    49896.51548529411764705882352941176470588
  136.  
  137. > p1(Pi);
  138.                                      446229
  139.                                      ------
  140.                                        11
  141.  
  142. > evalf(");
  143.                    40566.27272727272727272727272727272727273
  144.  
  145. > p2(Pi);
  146.                                     7689337
  147.                                     -------
  148.                                        85
  149.  
  150. > evalf(");
  151.                    90462.78823529411764705882352941176470588
  152.  
  153. > p(Pi/2);
  154.                / 223         1/2                   1/2      1/2\2
  155.                |----- I 44000    + 1/680 3392963053    68000   |
  156.                \44000                                          /
  157.  
  158. # This is about as small as the function gets
  159. > evalc(evalf("));
  160.                                                            9
  161.                .4989651537227366310160427807486631016040*10
  162.  
  163.                     + 47494.53383446247774655548050782504198492 I
  164.  
  165. > p(Pi);
  166.                  /223       1/2   1/2    20         1/2   1/2\2
  167.                  |--- 446229    11    + ---- 7689337    85   |
  168.                  \ 11                    17                  /
  169.  
  170. > evalf(");
  171.                                                              10
  172.                  .5623743076386678984956221653763304287514*10
  173.  
  174. > quit
  175.  
  176. Kate (with thanks to Dave Hare and Greg Fee)
  177. Waterloo Maple Software Technical Support
  178. support@maplesoft.on.ca
  179.