Re: Trouble with coupled quadratic equations where the solutions are degenerate/symmetric
- To: mathgroup at smc.vnet.net
- Subject: [mg108256] Re: [mg106476] Trouble with coupled quadratic equations where the solutions are degenerate/symmetric
- From: Robert Hoy <robert.hoy at yale.edu>
- Date: Fri, 12 Mar 2010 07:07:48 -0500 (EST)
- References: <201001141046.FAA19729@smc.vnet.net> <4B4F4A54.1060707@wolfram.com>
Hello - sorry for the extremely belated reply. Daniel, your method
worked, of course - thanks.
What I ended up going with was a numerical solver that makes initial
guesses (with some randomnness) as to the values of the {x, y, z}, and
then uses Newton iterations to try and solve them. If it doesn't work
after a number of iterations (which I set), it makes another initial
guess. The code uses LU decomposition and a number of standard
tricks. One thing I thought worth mentioning is, the method does not
employ complex numbers, and I don't see why complex numbers would be
needed. So, bringing things back to Mathematica, three questions:
1) Is there a way to restrict the Mathematica solver to work only in
real space?
2) A problem with solving these sets of equations in Mathematica is
its use of Grobner bases; solutions slow drastically as the number of
variables (3N) exceeds 21 (N = 7). Does Mathematica have a special
solver for quadratic forms (my equations are all quadratic forms)
which uses something other than Grobner bases?
3) To be clear, though I am solving equations in 3N variables, the
solutions 'live in' 3-dimensional Cartesian space (R^3). Is there
some command to tell Mathematica to look for solutions only in R^3?
The reason I ask is, the numerical solver I'm using also starts to
break as N becomes 'large' (how large 'large' is can be increased by
using tricks, but only up to a point), but I guess Mathematica might
include some special routines for dealing with this problem (which is,
after all, just solving real quadratic forms).
Thanks,
Robert Hoy
On Jan 14, 2010, at 11:46 AM, Daniel Lichtblau wrote:
> Robert Hoy wrote:
>> Hi, all! I am trying to solve for the structures of small atomic
>> clusters using Mathematica. Using adjacency (contact) matrices as
>> input, my the desired output is 3-D Cartesian coordinates. I
>> formulate the problem as coupled quadratic equations and
>> inequalities. For example, for clusters of M = 6 atoms I
>> introduce the variables
>> xvec = {{x1, y1, z1}, {x2, y2, z2}, {x3, y3, z3}, {x4, y4, z4},
>> {x5, y5,
>> z5}, {x6, y6, z6}}
>> Then, for an adjacency matrix A
>> {
>> {0, 1, 0, 1, 1, 1},
>> {1, 0, 1, 1, 0, 1},
>> {0, 1, 0, 1, 1, 1},
>> {1, 1, 1, 0, 1, 0},
>> {1, 0, 1, 1, 0, 1},
>> {1, 1, 1, 0, 1, 0}
>> }
>> there are M(M-1)/2 coupled equations and inequalities:
>> EqnTab = {(x1 - x2)^2 + (y1 - y2)^2 + (z1 - z2)^2 ==
>> 1, (x1 - x3)^2 + (y1 - y3)^2 + (z1 - z3)^2 >
>> 1, (x1 - x4)^2 + (y1 - y4)^2 + (z1 - z4)^2 ==
>> 1, (x1 - x5)^2 + (y1 - y5)^2 + (z1 - z5)^2 ==
>> 1, (x1 - x6)^2 + (y1 - y6)^2 + (z1 - z6)^2 ==
>> 1, (x2 - x3)^2 + (y2 - y3)^2 + (z2 - z3)^2 ==
>> 1, (x2 - x4)^2 + (y2 - y4)^2 + (z2 - z4)^2 ==
>> 1, (x2 - x5)^2 + (y2 - y5)^2 + (z2 - z5)^2 >
>> 1, (x2 - x6)^2 + (y2 - y6)^2 + (z2 - z6)^2 ==
>> 1, (x3 - x4)^2 + (y3 - y4)^2 + (z3 - z4)^2 ==
>> 1, (x3 - x5)^2 + (y3 - y5)^2 + (z3 - z5)^2 ==
>> 1, (x3 - x6)^2 + (y3 - y6)^2 + (z3 - z6)^2 ==
>> 1, (x4 - x5)^2 + (y4 - y5)^2 + (z4 - z5)^2 ==
>> 1, (x4 - x6)^2 + (y4 - y6)^2 + (z4 - z6)^2 >
>> 1, (x5 - x6)^2 + (y5 - y6)^2 + (z5 - z6)^2 == 1}
>> In these equations and inequalities, the left hand sides are
>> Dot[xvec[[m]] - xvec[[n]], xvec[[m]] - xvec[[n]] (i. e. the
>> squared Euclidean distance between xvec[[m]] and xvec[[n]] and the
>> 'right hand sides' are '==1' if A[[m,n]]==1 and '> 1' if
>> A[[m,n]]=0. Only terms above the diagonal of the adjacency matrix
>> (n > m) are necessary since A is symmetric.
>> I know the solution to these equations; the six elements of xvec
>> form the vertices of an octahedron.
>> Mathematica hangs when I try to solve this set. This isn't
>> surprising, because the above equations alone are invariant under
>> translations and rotations, so the solutions are continuously
>> degenerate. To eliminate the continuous degeneracies I add
>> equations constraining the first three 'atoms' to a plane (I can
>> do this because they are noncollinear). There is another discrete
>> mirror degeneracy (the other atoms can lie on either side of the
>> defined plane), so I add another inequality saying the fourth atom
>> is on one side. This adds ten equations/inequalities and gives
>> the system
>> CompleteEqnTab = {x1 == 0, y1 == 0, z1 > 0, x2 == 0, z2 == 0, y2 >
>> 0, y3 == 0, z3 == 0, x3 > 0,
>> x4/x3 + y4/y2 + z4/z1 > 1, [the 15 terms of EqnTab]}
>> Now, the problem is, I have been unable (after a considerable
>> amount of effort) to solve this system with Mathematica.
>> Reduce[CompleteEqnTab, Flatten[xvec], Reals] fails. So does
>> FindInstance[CompleteEqnTab, Flatten[xvec], Reals]. The memory
>> used by the kernel just increases (well past 1GB) until the
>> program hangs.
>> I think the problem is that the solution for xvec (i. e. the
>> vertices of an octahedron) has many rotational symmetries - its
>> point group is O_h. Perhaps these symmetries give Mathematica's
>> solver trouble?
>> In contrast, the exact same approach for a less symmetric cluster
>> with adjacency matrix
>> (
>> {0, 1, 1, 1, 1, 1},
>> {1, 0, 1, 1, 1, 1},
>> {1, 1, 0, 1, 0, 0},
>> {1, 1, 1, 0, 1, 0},
>> {1, 1, 0, 1, 0, 1},
>> {1, 1, 0, 0, 1, 0}
>> )
>> and equation/inequality set
>> CompleteEqnTab2 = {x1 == 0, y1 == 0, z1 > 0, x2 == 0, z2 == 0, y2 >
>> 0, y3 == 0, z3 == 0, x3 > 0,
>> x4/x3 + y4/y2 + z4/z1 > 1, (x1 - x2)^2 + (y1 - y2)^2 + (z1 -
>> z2)^2 ==
>> 1, (x1 - x3)^2 + (y1 - y3)^2 + (z1 - z3)^2 ==
>> 1, (x1 - x4)^2 + (y1 - y4)^2 + (z1 - z4)^2 ==
>> 1, (x1 - x5)^2 + (y1 - y5)^2 + (z1 - z5)^2 ==
>> 1, (x1 - x6)^2 + (y1 - y6)^2 + (z1 - z6)^2 ==
>> 1, (x2 - x3)^2 + (y2 - y3)^2 + (z2 - z3)^2 ==
>> 1, (x2 - x4)^2 + (y2 - y4)^2 + (z2 - z4)^2 ==
>> 1, (x2 - x5)^2 + (y2 - y5)^2 + (z2 - z5)^2 ==
>> 1, (x2 - x6)^2 + (y2 - y6)^2 + (z2 - z6)^2 ==
>> 1, (x3 - x4)^2 + (y3 - y4)^2 + (z3 - z4)^2 ==
>> 1, (x3 - x5)^2 + (y3 - y5)^2 + (z3 - z5)^2 >
>> 1, (x3 - x6)^2 + (y3 - y6)^2 + (z3 - z6)^2 >
>> 1, (x4 - x5)^2 + (y4 - y5)^2 + (z4 - z5)^2 ==
>> 1, (x4 - x6)^2 + (y4 - y6)^2 + (z4 - z6)^2 >
>> 1, (x5 - x6)^2 + (y5 - y6)^2 + (z5 - z6)^2 == 1}
>> gives a quick solution
>> In[32]:= Timing[Reduce[CompleteEqnTab2, Flatten[xvec], Reals]]
>> Out[32]= {18.527, x1 == 0 && y1 == 0 && z1 == 1/Sqrt[2] && x2 == 0
>> && y2 == 1/Sqrt[2] &&
>> z2 == 0 && x3 == 1/Sqrt[2] && y3 == 0 && z3 == 0 && x4 == 1/
>> Sqrt[2] &&
>> y4 == 1/Sqrt[2] && z4 == 1/Sqrt[2] && x5 == -(1/(3 Sqrt[2])) &&
>> y5 == (2 Sqrt[2])/3 && z5 == (2 Sqrt[2])/3 && x6 == -(11/(9
>> Sqrt[2])) &&
>> y6 == 5/(9 Sqrt[2]) && z6 == 5/(9 Sqrt[2])}
>> Similar issues hold for larger clusters; Mathematica solves the
>> less- symmetric systems and hangs on the more-symmetric ones. Is
>> there anything I can do about this? I thought FindInstance was
>> designed for this sort of thing but it doesn't seem to help.
>> As you can see from the above, the form of the equations and
>> inequalities is simple. I assume Mathematica can handle the
>> symmetric cases and I'm just unaware of the necessary
>> functionality. Any help would be greatly appreciated. I am using
>> Mathematica 7.0.1.0 on Mac OS 10.5.8.
>> Thanks,
>> Robert Hoy
>
> I'm not sure what I am trying is quite what you need. Might give a
> start though. The idea is this.
>
> (1) Start with the equations (that is, ignore the inequalities)
> (2) Reduce by the variables that are set to zero,
> (3) Perturb numerically
> (4) Solve numerically
> (5) Discard complex solutions
> (6) Put back the zero-valued variables
> (7) Apply the inequalities
> (8) Inspect what remains as solution, if anything. If it gives some
> insight into the exact solutions, use that to further simplify
> matters.
>
> I show this below.
>
> In[53]:= CompleteEqnTab =
> Join[{x1 == 0, y1 == 0, z1 > 0, x2 == 0, z2 == 0, y2 > 0, y3 == 0,
> z3 == 0, x3 > 0, x4/x3 + y4/y2 + z4/z1 > 1}, EqnTab];
>
> In[54]:= eqns = Cases[CompleteEqnTab, _Equal];
>
> In[55]:= polys = Apply[Subtract, eqns, 1];
>
> In[56]:= rules = Thread[{x1, y1, x2, z2, y3, z3} -> 0];
>
> In[57]:= p2 = DeleteCases[polys /. rules, 0]
>
> Out[57]= {-1 + y2^2 + z1^2, -1 + x4^2 + y4^2 + (z1 - z4)^2, -1 +
> x5^2 + y5^2 + (z1 - z5)^2, -1 + x6^2 + y6^2 + (z1 - z6)^2, -1 +
> x3^2 + y2^2, -1 + x4^2 + (y2 - y4)^2 + z4^2, -1 +
> x6^2 + (y2 - y6)^2 + z6^2, -1 + (x3 - x4)^2 + y4^2 +
> z4^2, -1 + (x3 - x5)^2 + y5^2 + z5^2, -1 + (x3 - x6)^2 + y6^2 +
> z6^2, -1 + (x4 - x5)^2 + (y4 - y5)^2 + (z4 - z5)^2, -1 + (x5 -
> x6)^2 + (y5 - y6)^2 + (z5 - z6)^2}
>
> In[58]:= vars2 = Variables[p2];
>
> In[59]:= p3 = p2 + 1/10000*RandomReal[{0, 1}, Length[p2]]
>
> Out[59]= {-0.999914 + y2^2 + z1^2, -0.999904 + x4^2 +
> y4^2 + (z1 - z4)^2, -0.999945 + x5^2 +
> y5^2 + (z1 - z5)^2, -0.999912 + x6^2 +
> y6^2 + (z1 - z6)^2, -0.999905 + x3^2 + y2^2, -0.999958 +
> x4^2 + (y2 - y4)^2 + z4^2, -0.999909 + x6^2 + (y2 - y6)^2 +
> z6^2, -0.999983 + (x3 - x4)^2 + y4^2 +
> z4^2, -0.999903 + (x3 - x5)^2 + y5^2 +
> z5^2, -0.999948 + (x3 - x6)^2 + y6^2 +
> z6^2, -0.999968 + (x4 - x5)^2 + (y4 - y5)^2 + (z4 -
> z5)^2, -0.999984 + (x5 - x6)^2 + (y5 - y6)^2 + (z5 - z6)^2}
>
> In[60]:= Timing[solns = NSolve[p3];]
>
> Out[60]= {29.8675, Null}
>
> In[61]:= s2 = vars2 /. solns;
>
> In[62]:= keep = Map[Apply[And, Map[Head[#] =!= Complex &, #]] &, s2];
>
> In[63]:= realsols = Pick[solns, keep, True];
>
> In[64]:= realsols = Map[Join[rules, #] &, realsols];
>
> ineqs = DeleteCases[CompleteEqnTab, _Equal];
>
> In[66]:= keep2 = Apply[And, ineqs /. realsols, 1]
>
> Out[66]= {False, True, False, False, False, False, False, False, \
> False, False, False, False, False, False, False, False, False,
> False, \
> False, False, False, False, False, False, False, False, False,
> False, \
> False, False, False, False, False, False, False, False, False,
> False, \
> False, False, False, False, False, False, False, False, False, False}
>
> In[67]:= Pick[realsols, keep2, True]
>
> Out[67]= {{x1 -> 0, y1 -> 0, x2 -> 0, z2 -> 0, y3 -> 0, z3 -> 0,
> y2 -> 0.00324073, x4 -> 0.502247, x6 -> 0.497657, z4 -> 0.502289,
> y5 -> -0.00324273, z6 -> 0.497678, y4 -> 0.707095, y6 -> -0.707086,
> z1 -> 0.999952, z5 -> 0.999946, x3 -> 0.999947, x5 -> 0.999967}}
>
> That finishes steps (1)-(7). For (8), I would surmise that you want
> to impose the following equations: {z1==1, z5==1, x3==1, x5==1,
> x4==z4, y4==-y6, x6==z6, y2==-y5}.
>
> In[75]:= moreeqns = {z1 == 1, z5 == 1, x3 == 1, x5 == 1, x4 == z4,
> y4 == -y6, x6 == z6, y2 == -y5};
>
> In[80]:= fulleqns = Join[eqns, moreeqns];
>
> In[77]:= fullpolys = Apply[Subtract, fulleqns, 1];
>
> In[78]:= Timing[solns2 = NSolve[fullpolys];]
>
> Out[78]= {0.110983, Null}
>
> Let's have a look.
>
> In[79]:= solns2
>
> Out[79]= {{x1 -> 0., x3 -> 1., y1 -> 0., y3 -> 0., z1 -> 1., z3 -> 0.,
> z2 -> 0., z5 -> 1., z4 -> 0.5, z6 -> 0.5, y4 -> -0.707107,
> y6 -> 0.707107, y2 -> 0., y5 -> 0., x4 -> 0.5, x6 -> 0.5, x2 -> 0.,
> x5 -> 1.}, {x1 -> 0., x3 -> 1., y1 -> 0., y3 -> 0., z1 -> 1.,
> z3 -> 0., z2 -> 0., z5 -> 1., z4 -> 0.5, z6 -> 0.5, y4 -> -0.707107,
> y6 -> 0.707107, y2 -> 0., y5 -> 0., x4 -> 0.5, x6 -> 0.5, x2 -> 0.,
> x5 -> 1.}, {x1 -> 0., x3 -> 1., y1 -> 0., y3 -> 0., z1 -> 1.,
> z3 -> 0., z2 -> 0., z5 -> 1., z4 -> 0.5, z6 -> 0.5, y4 -> 0.707107,
> y6 -> -0.707107, y2 -> 0., y5 -> 0., x4 -> 0.5, x6 -> 0.5, x2 -> 0.,
> x5 -> 1.}, {x1 -> 0., x3 -> 1., y1 -> 0., y3 -> 0., z1 -> 1.,
> z3 -> 0., z2 -> 0., z5 -> 1., z4 -> 0.5, z6 -> 0.5, y4 -> 0.707107,
> y6 -> -0.707107, y2 -> 0., y5 -> 0., x4 -> 0.5, x6 -> 0.5, x2 -> 0.,
> x5 -> 1.}}
>
> Not hard to deduce the algebraic values (all are in the set +-
> {0,1/2,Sqrt[2]/2,1}). But they all seem to violate the inequality
> that y2>0.
>
> Notice that the enforcement of real solutions might also be overly
> restrictive. Since we do a perturbation, there might be complex
> solutions with small imaginary parts that would be zero in the
> unperturbed case. So one might want to allow them, modifying some of
> the subsequent code as needed to check inequalities, etc. When I
> tried that I still seem to run afoul of the y2>0 requirement. You
> might have better luck.
>
> Daniel Lichtblau
> Wolfram Research