From: jcao@option.GBA.NYU.EDU (Jingbin Cao)
Newsgroups: sci.math.symbolic
Subject: Mma Solve[] oddity (a much shorter version)
Summary: Solve[] finds roots in a strang way, and may lose roots
Keywords: Solve, oddity
Date: 21 Aug 92 08:53:21 GMT
Organization: NYU Stern School of Business

The following log documents some weired Solve[] oddities. They are:

	1. Solve[] cannot find roots in some cases, at least not within
	   reasonable period of time. But if you divide the equations
	   by something, Solve[] finds some roots!
	2. Solve[] may lose some roots in the process of finding them,
	   even in some trivial cases. So there is no guarrenty that 
	   you can get all the roots.

I apologize if the log is too long, but I wanted to capture the whole
log.  Viewing it using some editor like emacs is recommended. My
comments quoted by (* *) were added after or during the log.

------------------------------log begin---------------------------------
Mathematica 2.0 for SPARC
Copyright 1988-91 Wolfram Research, Inc.
 -- Terminal graphics initialized -- 

In[1]:= u1=-b1 + b2 + 3*n - 3*n*q + 2*x1 + x1^2 - 2*x2 - x2^2

                                          2            2
Out[1]= -b1 + b2 + 3 n - 3 n q + 2 x1 + x1  - 2 x2 - x2

In[2]:= u2=-b1 + b2 - n + n*q - 2*x1 - 4*n*x1 + 4*n*q*x1 - 3*x1^2 + 2*x2 + 
       4*x1*x2 - x2^2

                                                            2
Out[2]= -b1 + b2 - n + n q - 2 x1 - 4 n x1 + 4 n q x1 - 3 x1  + 2 x2 + 
 
                 2
>    4 x1 x2 - x2

In[3]:= v1=b1 - b2 + 3*n - 3*n*q + 4*x1 - x1^2 - 4*x2 + x2^2

                                         2            2
Out[3]= b1 - b2 + 3 n - 3 n q + 4 x1 - x1  - 4 x2 + x2

In[4]:= v2=-b1 + b2 + 5*n - 5*n*q + 4*x1 + x1^2 - 4*x2 - 4*n*x2 + 4*n*q*x2 - 
       4*x1*x2 + 3*x2^2

                                          2
Out[4]= -b1 + b2 + 5 n - 5 n q + 4 x1 + x1  - 4 x2 - 4 n x2 + 4 n q x2 - 
 
                   2
>    4 x1 x2 + 3 x2

In[5]:= f=u1 u2;

In[6]:= g=v1 v2;

In[7]:= Timing[sol=Solve[{f==0, g==0}, {x1, x2}]]

	(* This is basically what I want to do: to find out all the
	   roots for {f==0, g==}. However, some weired things are revealed *)

Interrupt> a  (* It took more than a minute of cpu time,
		but the computation is not done. So I aborted it *)

Out[7]= $Aborted

In[8]:= w=18*(-n + n*q - x1 + x2)^2

                               2
Out[8]= 18 (-n + n q - x1 + x2)

In[9]:= f=f/w;

In[10]:= g=g/w;

In[11]:= Timing[sol=Solve[{f==0, g==0}, {x1, x2}]]

Out[11]= {19.4667 Second, {{x1 -> 
 
        -8 + 4 Sqrt[9 + 4 b1 - 4 b2 - 12 n + 12 n q]        1
>       --------------------------------------------, x2 -> -}, 
                             8                              2
 
             -8 - 4 Sqrt[9 + 4 b1 - 4 b2 - 12 n + 12 n q]        1
>     {x1 -> --------------------------------------------, x2 -> -}, 
                                  8                              2
 
             -(-15 - 4 b1 + 4 b2 + 20 n - 20 n q)   3 (12 - 16 n + 16 n q)
>     {x1 -> ------------------------------------ - ----------------------, 
                     4 (3 - 4 n + 4 n q)             8 (3 - 4 n + 4 n q)
 
             -(-15 - 4 b1 + 4 b2 + 20 n - 20 n q)
>      x2 -> ------------------------------------}, 
                     4 (3 - 4 n + 4 n q)
 
             1        16 + 4 Sqrt[9 - 4 b1 + 4 b2 - 12 n + 12 n q]
>     {x1 -> -, x2 -> --------------------------------------------}, 
             2                             8
 
             1        16 - 4 Sqrt[9 - 4 b1 + 4 b2 - 12 n + 12 n q]
>     {x1 -> -, x2 -> --------------------------------------------}}}
             2                             8

	(* So a little trick brings us five pairs of roots, and it only
	   took ~19.5 cpu seconds.
	   Why dividing f and g by w helps is beyond my comprehension.
	   Now let's have a closer look. *)

In[12]:= Table[Factor[u1/.sol[[i]]], {i, 5}]

                (-9 - 4 b1 + 4 b2 + 12 n - 12 n q) (3 - 2 n + 2 n q)
Out[12]= {0, 0, ----------------------------------------------------, 
                                2 (3 - 4 n + 4 n q)
 
     -36 + 24 n - 24 n q - 12 Sqrt[9 - 4 b1 + 4 b2 - 12 n + 12 n q]
>    --------------------------------------------------------------, 
                                   4
 
     -36 + 24 n - 24 n q + 12 Sqrt[9 - 4 b1 + 4 b2 - 12 n + 12 n q]
>    --------------------------------------------------------------}
                                   4

In[13]:= Table[Factor[u2/.sol[[i]]], {i, 5}]

Out[13]= {(-9 - 4 b1 + 4 b2 + 12 n - 12 n q + 
 
>       12 Sqrt[9 + 4 b1 - 4 b2 - 12 n + 12 n q] - 
 
>       8 n Sqrt[9 + 4 b1 - 4 b2 - 12 n + 12 n q] + 
 
>       8 n q Sqrt[9 + 4 b1 - 4 b2 - 12 n + 12 n q] - 
 
>       3 (9 + 4 b1 - 4 b2 - 12 n + 12 n q)) / 4, 
 
>    (-9 - 4 b1 + 4 b2 + 12 n - 12 n q - 
 
>       12 Sqrt[9 + 4 b1 - 4 b2 - 12 n + 12 n q] + 
 
>       8 n Sqrt[9 + 4 b1 - 4 b2 - 12 n + 12 n q] - 
 
>       8 n q Sqrt[9 + 4 b1 - 4 b2 - 12 n + 12 n q] - 
 
>       3 (9 + 4 b1 - 4 b2 - 12 n + 12 n q)) / 4, 0, 0, 0}

In[14]:= Table[Factor[v1/.sol[[i]]], {i, 5}]

          -36 + 24 n - 24 n q + 12 Sqrt[9 + 4 b1 - 4 b2 - 12 n + 12 n q]
Out[14]= {--------------------------------------------------------------, 
                                        4
 
     -36 + 24 n - 24 n q - 12 Sqrt[9 + 4 b1 - 4 b2 - 12 n + 12 n q]
>    --------------------------------------------------------------, 
                                   4
 
     (3 - 2 n + 2 n q) (9 - 4 b1 + 4 b2 - 12 n + 12 n q)
>    ---------------------------------------------------, 0, 0}
                    2 (-3 + 4 n - 4 n q)

In[15]:= Table[Factor[v2/.sol[[i]]], {i, 5}]

Out[15]= {0, 0, 0, (9 - 4 b1 + 4 b2 - 12 n + 12 n q + 
 
>       12 Sqrt[9 - 4 b1 + 4 b2 - 12 n + 12 n q] - 
 
>       8 n Sqrt[9 - 4 b1 + 4 b2 - 12 n + 12 n q] + 
 
>       8 n q Sqrt[9 - 4 b1 + 4 b2 - 12 n + 12 n q] + 
 
>       3 (9 - 4 b1 + 4 b2 - 12 n + 12 n q)) / 4, 
 
>    (9 - 4 b1 + 4 b2 - 12 n + 12 n q - 
 
>       12 Sqrt[9 - 4 b1 + 4 b2 - 12 n + 12 n q] + 
 
>       8 n Sqrt[9 - 4 b1 + 4 b2 - 12 n + 12 n q] - 
 
>       8 n q Sqrt[9 - 4 b1 + 4 b2 - 12 n + 12 n q] + 
 
>       3 (9 - 4 b1 + 4 b2 - 12 n + 12 n q)) / 4}

(* So,
	sol[[1]], sol[[2]] are solutions to {u1==0, v2==0},
	sol[[3]] 	   is  solution  to {u2==0, v2==0}, and
	sol[[4]], sol[[5]] are solutions to {u2==0, v1==0}.

   Now, some more questions:
	Q1: Is there solution to {u1==0, v1==0}?
	Q2: Is there _other_ solution to {u2==0, v2==0}?
	Q3: Is there _other_ solution to {u1==0, v2==0} or {u2==0, v1==0}?
*)

In[16]:= Timing[sol1=Solve[{u1==0, v1==0}, {x1, x2}]] (* for Q1 *)

Out[16]= {0.4 Second, {{x1 -> 
 
                                                      2            2      2  2
        -((n - n q) (-2 n + 2 n q))   -b1 + b2 + n + n  - n q - 2 n  q + n  q
>       --------------------------- - ----------------------------------------
               2 (-n + n q)                         2 (-n + n q)
 
                                  2            2      2  2
                -(-b1 + b2 + n + n  - n q - 2 n  q + n  q )
>       , x2 -> -------------------------------------------}}}
                               2 (-n + n q)

	(* So we find one more pair of roots for {u1==0, v1==0} *)

In[17]:= Timing[sol1=Solve[{u1==0, v2==0}, {x1, x2}]] (* for Q2 *)

Interrupt> a  (* I have to abort it again because it had taken
		more than 7 cpu minutes *)

Out[17]= $Aborted

	(* Now let's look at Q3 *)

In[18]:= Timing[sol1=Solve[{u2==0, v2==0}, {x1, x2}]] 

Out[18]= {31.8333 Second, {{x1 -> 
 
                                          9
>       (-5 + 10 n - 10 n q + ------------------------- - 
                              (1 - q) (3 - 4 n + 4 n q)
 
	......................

	(* the output is way too long and I have to delete it.
	   the roots are simplified in In[21] *)

In[19]:= sol2=sol1; sol1=%16; (* let me correct my mistake first *)

In[20]:= Table[{Factor[x1/.sol2[[i]]], factor[x2/.sol2[[i]]]}, {i, 2}]

General::spell1: 
   Possible spelling error: new symbol name "factor"
     is similar to existing symbol "Factor".

Interrupt> a     (* opps, another mistake *)

Out[20]= $Aborted

In[21]:= Table[{Factor[x1/.sol2[[i]]], Factor[x2/.sol2[[i]]]}, {i, 2}]

           -3 + 4 b1 - 4 b2 + 4 n - 4 n q
Out[21]= {{------------------------------, 
                4 (3 - 4 n + 4 n q)
 
      -15 - 4 b1 + 4 b2 + 20 n - 20 n q
>     ---------------------------------}, 
            4 (-3 + 4 n - 4 n q)
 
                     2            2      2  2
      b1 - b2 - n + n  + n q - 2 n  q + n  q
>    {---------------------------------------, 
                   2 n (-1 + q)
 
                      2            2      2  2
      -b1 + b2 + n + n  - n q - 2 n  q + n  q
>     ----------------------------------------}}
                    2 n (1 - q)

In[22]:= (* So we find one more pair of roots for {u2==0, v2==0}

	What are my conclusions of this exercise?

	1. Solve[] cannot find roots in some cases, at least not within
	   reasonable period of time. But if you divide the equations
	   by something, Solve[] finds some roots!
	2. Solve[] may lose some roots in the process of finding them,
	   even in some trivial cases. *)
---------------------------------log ends---------------------------------
Newsgroups: sci.math.symbolic