Re: Trouble with coupled quadratic equations where the
- To: mathgroup at smc.vnet.net
- Subject: [mg106512] Re: [mg106476] Trouble with coupled quadratic equations where the
- From: Daniel Lichtblau <danl at wolfram.com>
- Date: Fri, 15 Jan 2010 03:17:34 -0500 (EST)
- References: <201001141046.FAA19729@smc.vnet.net>
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
- References:
- Trouble with coupled quadratic equations where the solutions are degenerate/symmetric
- From: Robert Hoy <robert.hoy@yale.edu>
- Trouble with coupled quadratic equations where the solutions are degenerate/symmetric