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