Re: Trouble with coupled quadratic equations where the solutions are
- To: mathgroup at smc.vnet.net
- Subject: [mg106617] Re: Trouble with coupled quadratic equations where the solutions are
- From: Ray Koopman <koopman at sfu.ca>
- Date: Mon, 18 Jan 2010 05:40:01 -0500 (EST)
- References: <himsmg$j8f$1@smc.vnet.net>
On Jan 14, 2:46 am, Robert Hoy <robert.... at yale.edu> 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.
The squared distances corresponding to the off-diagonal zeros in A
should all be 2. If we change the "> 1"s in EqnTab to "== 2"s, then
Reduce should be able to solve it, but it doesn't. Instead, it returns
"False" almost immediately. However, if we swap the labels of atoms 3
and 4, so that A becomes
{{0, 1, 1, 0, 1, 1},
{1, 0, 1, 1, 0, 1},
{1, 1, 0, 1, 1, 0},
{0, 1, 1, 0, 1, 1},
{1, 0, 1, 1, 0, 1},
{1, 1, 0, 1, 1, 0}},
then Reduce solves it both when the "> 1"s are changed to "== 2"s,
and when the first "> 1" is left as is but the other two are changed
to "== (x1 - x4)^2 + (y1 - y4)^2 + (z1 - z4)^2" to force the three
distances to be equal. However, Reduce still fails to solve it when
all three distances are simply "> 1".
>
> 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