From: rusin@vesuvius.math.niu.edu (Dave Rusin) Newsgroups: sci.math.symbolic Subject: Re: Help! Date: 10 Mar 1998 20:03:23 GMT In article <3504D5D9.54B6BFB1@adam.montreal.ca>, root wrote: >I am trying to find solutions to the following system of equations. > >Let g[u1, u2, q]= u1^8 (u1 - q^2 u2)(u1 u2 - q^2 ) - (u2 - q^2 u2)(1 - >q^2 u1 u2) > >The system is then >{g[u1, u2, q] = 0, g[u2, u1, q] =0} Let me analyze the solution in some detail, in case I've misunderstood what your real goal is. Your two equations describe a curve in u1-u2-q space. I'm not sure what exactly you mean by "find the solutions", but I'm guessing you want to know u1 and u2 as functions of q. Your other comments suggest you want to be able to compute them numerically for given numerical q. (If that's true, I can't see why you refer to q since only Q:= q^2 is involved. I will refer only to Q in the sequel.) Maybe it's easier if you first look at the projection of your curve to u1-u2 space: I used maple. >readlib(eliminate): >eliminate({g,h},Q); This gives a formula expression Q as a rational function of u1, u2, and an equation describing the curve in the u1-u2 plane. That curve includes a few easy components (u1= +- u2 and u1 = +- 1/u2) as well as one complicated one: 14 14 7 15 15 7 11 11 9 13 13 9 9 15 u1 + u2 + u2 u1 + u2 u1 + u2 u1 + u2 u1 + u2 u1 - u2 u1 11 13 9 15 11 13 13 2 14 6 10 - u2 u1 - u1 u2 - u1 u2 - u2 u1 - 2 u1 u2 - 4 u2 u1 12 4 8 8 10 6 12 4 11 3 - 3 u1 u2 - 5 u2 u1 - 4 u2 u1 - 3 u2 u1 - 2 u1 u2 6 8 8 6 10 4 15 9 5 2 12 + 4 u2 u1 + 4 u2 u1 + 3 u2 u1 + u1 u2 - 3 u2 u1 + 2 u1 u2 4 10 12 2 11 5 3 13 9 7 + 3 u2 u1 + 2 u1 u2 + 3 u1 u2 + 2 u1 u2 + 4 u2 u1 7 9 3 11 5 11 16 16 13 15 + 4 u2 u1 - 2 u1 u2 + 3 u1 u2 - u2 - u1 - u1 u2 + u2 u1 13 3 14 2 7 7 5 9 + 2 u1 u2 - 2 u1 u2 - 4 u2 u1 - 3 u2 u1 = 0 (The equation Q=... then allows one to reconstruct the curve in 3-dimensional space as a covering of this curve in the plane.) I'll discuss this complicated curve carefully; however, we will eventually discover it contains _none of the points of interest in the original problem_. I observe that maple's implicitplot feature has some difficulty discerning the behaviour of this curve in a few regions; let's analyze the curve carefully. In order to visualize this planar curve, we can convert to polar coordinates. I will avoid trigonometry for a moment: > subs({u1=t*v1,u2=t*v2}, crv): > simplify("/t^14): > simplify(", {t^2=T}): > c2:=convert(factor(taylor(",T,6)),polynom); which gives a quintic in T: 2 2 2 2 6 3 3 6 4 4 c2 := (v2 - v2 v1 + v1 ) (v2 + v1 ) (v2 - v1 v2 + v1 ) (v2 + v1 ) - ... 9 9 2 2 4 4 5 - v2 v1 (v2 + v1 ) (v2 + v1 ) T Thus it shouldn't be hard to draw this curve in polar coordinates: for any point (v1,v2) on the unit circle, there should be at most 5 points on the curve which lie on the ray through that point. Actually the number of points per ray is locally constant since there is no bifurcation of the roots: > c3:=subs({v1=cos(H),v2=sin(H)},c2): > eliminate({c3,diff(c3,T)},T); yields an immense polynomial in cos(H), sin(H) which never vanishes. On the other hand, the polynomial drops degree precisely at the axes, so the number of roots can change there. Taking a few test values of H shows there is one point per ray in the first and third quadrants and two per ray in the other quadrants. (There are three real roots T in those quadrants, but only two are positive, so t is one of two pairs of positive/negative values.) Thus we can ask maple to plot some points, e.g.: > for ii from 1 to 100 do fsolve(evalf(subs(H=ii*Pi/400,c3))); od; (By symmetry of the curve we need only the T for H in [-Pi/4, Pi/4]). The values of T rise monotonically from 1 to 2 in this interval. As H varies from 0 to -Pi/4, one values of T rises from 1 to about 1.568670203; the other positive one decreases (from +infinity) to about 4.096085035. A plot of these points shows a closed loop which roughly resembles the outline of a rounded square with one pair of opposite corners pulled away from the square. The unbounded portion of the curve resembles a hyperbola xy=-2. The bounded and unbounded portions of this curve describe the "nontrivial" combinations of u1 and u2 for which a Q exists satisfying the original equations g=h=0. As I remarked at the beginning, there are also the curves u1=u2, u1=-u2, u1=1/u2, and u1=-1/u2. Now, this is just the projection of the original curve to the u1-u2 plane. The original curve in R^3 is easily found from this since we know how to lift the projection back up. Indeed, elimination of Q at the beginning provided a rational function which enables us to write Q = 4 4 4 3 2 2 3 4 (v2 + v1 ) - v2 v1 (v2 - v2 v1 + v2 v1 - v2 v1 + v1 ) 4 3 2 2 3 4 2 8 8 5 (v2 + v2 v1 + v2 v1 + v2 v1 + v1 ) T + v2 v1 T ---------------------------------------------------------------- 4 4 2 2 6 3 3 6 (v2 + v1 ) - (v2 - v2 v1 + v1 ) (v2 - v1 v2 + v1 ) T 2 2 4 4 2 7 7 4 - v2 v1 (v2 + v1 ) T + v2 v1 T where the dependence of T on the point (v1,v2) on the unit circle has been discussed above. Here the denominator is only zero when v1=v2 or v1 v2 = 0 (as can be seen by eliminating T between the denominator of Q and the equation c2 which defines T). In the former case, the numerator is also zero, and in fact the expression has a limit of 1; the curve in 3 dimensions is locally isomorphic to its projection in R^2 here. On the other hand, as either v_i -> 0 we get Q -> +- oo (e.g. the near H=0, Q is (1/H)(1 - H^2/3 + ...) We can similarly show Q is never zero. In particular, Q never changes sign along the curve except at those four asymptotes. A quick test of a few values shows that Q is negative on the unbounded portion of the curve, and on the bounded portion when u1 u2 < 0. In the original application, Q=q^2 is positive, we see (u1, u2) will have to lie in the first or third quadrants and on the bounded portion of the curve. However, in those regions, we find the minimum value of Q to be 1, so since -1 < q was originally assumed, we conclude * No points of interest in the original problem lie on the curve * discussed so far; all must instead project to the "trivial" * curves u1 = +- u2 or u1 = +- 1/u2. So now we return to the original problem and look for solutions with these extra conditions. If u1=1/u2, the only solutions are those with Q=1 (in which case u2 is arbitrary). If u1=-1/u2, the only solutions are those with Q=-1, in which case again u2 is arbitrary. If u1=-u2, then Q is now variable, but u2 must be determined from 2 8 10 (- 1 + Q) + Q (- 1 + Q) u2 + Q (1 + Q) u2 + (1 + Q) u2 = 0 If u1=u2, then either Q=1 (u2 arbitrary), u2=0, 1, or -1 (Q arbitrary), or 8 6 6 4 4 2 2 u2 + u2 - u2 Q + u2 - u2 Q + u2 - u2 Q + 1 = 0 Perhaps Mathematica was misbehaving because the algebraic variety in question has several distinct components, and it was unable to figure out how to pick the right one on which to perform its business. dave