Re: Trouble with coupled quadratic equations where the
Daniel Lichtblau
> 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
Artur
Who know how to execute SymbolicPolynomialMod
PolynomialMod working only on numbers (not for symbols)
If we know that a,b,c,d,f are five roots of general quintic polynomial
X^5+t X^4+p X^3+q X^2+r X+s==0
Then we can do SymbolicPolynomialMod of g (manually not by Mathemathica
unfortunatelly)
g=(a^6 b^7 c^5 d^2 + a^7 b^5 c^6 d^2 + a^5 b^6 c^7 d^2 +
a^7 b^6 c^2 d^5 + a^2 b^7 c^6 d^5 + a^6 b^2 c^7 d^5 +
a^5 b^7 c^2 d^6 + a^7 b^2 c^5 d^6 + a^2 b^5 c^7 d^6 +
a^6 b^5 c^2 d^7 + a^2 b^6 c^5 d^7 + a^5 b^2 c^6 d^7 -
a^6 b^8 c^3 d^2 f - a^8 b^3 c^6 d^2 f - a^3 b^6 c^8 d^2 f -
a^8 b^6 c^2 d^3 f - a^2 b^8 c^6 d^3 f - a^6 b^2 c^8 d^3 f -
a^3 b^8 c^2 d^6 f - a^8 b^2 c^3 d^6 f - a^2 b^3 c^8 d^6 f -
a^6 b^3 c^2 d^8 f - a^2 b^6 c^3 d^8 f - a^3 b^2 c^6 d^8 f +
a^7 b^6 c^5 f^2 + a^5 b^7 c^6 f^2 + a^6 b^5 c^7 f^2 -
a^8 b^6 c^3 d f^2 - a^3 b^8 c^6 d f^2 - a^6 b^3 c^8 d f^2 -
a^6 b^8 c d^3 f^2 - a^8 b c^6 d^3 f^2 - a b^6 c^8 d^3 f^2 +
a^6 b^7 d^5 f^2 + a^7 c^6 d^5 f^2 + b^6 c^7 d^5 f^2 +
a^7 b^5 d^6 f^2 - a^8 b^3 c d^6 f^2 - a b^8 c^3 d^6 f^2 +
b^7 c^5 d^6 f^2 + a^5 c^7 d^6 f^2 - a^3 b c^8 d^6 f^2 +
a^5 b^6 d^7 f^2 + a^6 c^5 d^7 f^2 + b^5 c^6 d^7 f^2 -
a^3 b^6 c d^8 f^2 - a^6 b c^3 d^8 f^2 - a b^3 c^6 d^8 f^2 -
a^6 b^8 c^2 d f^3 - a^8 b^2 c^6 d f^3 - a^2 b^6 c^8 d f^3 -
a^8 b^6 c d^2 f^3 - a b^8 c^6 d^2 f^3 - a^6 b c^8 d^2 f^3 -
a^2 b^8 c d^6 f^3 - a^8 b c^2 d^6 f^3 - a b^2 c^8 d^6 f^3 -
a^6 b^2 c d^8 f^3 - a b^6 c^2 d^8 f^3 - a^2 b c^6 d^8 f^3 +
a^6 b^7 c^2 f^5 + a^7 b^2 c^6 f^5 + a^2 b^6 c^7 f^5 +
a^7 b^6 d^2 f^5 + b^7 c^6 d^2 f^5 + a^6 c^7 d^2 f^5 +
a^2 b^7 d^6 f^5 + a^7 c^2 d^6 f^5 + b^2 c^7 d^6 f^5 +
a^6 b^2 d^7 f^5 + b^6 c^2 d^7 f^5 + a^2 c^6 d^7 f^5 +
a^7 b^5 c^2 f^6 + a^2 b^7 c^5 f^6 + a^5 b^2 c^7 f^6 -
a^8 b^3 c^2 d f^6 - a^2 b^8 c^3 d f^6 - a^3 b^2 c^8 d f^6 +
a^5 b^7 d^2 f^6 - a^3 b^8 c d^2 f^6 - a^8 b c^3 d^2 f^6 +
a^7 c^5 d^2 f^6 + b^5 c^7 d^2 f^6 - a b^3 c^8 d^2 f^6 -
a^8 b^2 c d^3 f^6 - a b^8 c^2 d^3 f^6 - a^2 b c^8 d^3 f^6 +
a^7 b^2 d^5 f^6 + b^7 c^2 d^5 f^6 + a^2 c^7 d^5 f^6 +
a^2 b^5 d^7 f^6 + a^5 c^2 d^7 f^6 + b^2 c^5 d^7 f^6 -
a^2 b^3 c d^8 f^6 - a^3 b c^2 d^8 f^6 - a b^2 c^3 d^8 f^6 +
a^5 b^6 c^2 f^7 + a^6 b^2 c^5 f^7 + a^2 b^5 c^6 f^7 +
a^6 b^5 d^2 f^7 + b^6 c^5 d^2 f^7 + a^5 c^6 d^2 f^7 +
a^2 b^6 d^5 f^7 + a^6 c^2 d^5 f^7 + b^2 c^6 d^5 f^7 +
a^5 b^2 d^6 f^7 + b^5 c^2 d^6 f^7 + a^2 c^5 d^6 f^7 -
a^3 b^6 c^2 d f^8 - a^6 b^2 c^3 d f^8 - a^2 b^3 c^6 d f^8 -
a^6 b^3 c d^2 f^8 - a b^6 c^3 d^2 f^8 - a^3 b c^6 d^2 f^8 -
a^2 b^6 c d^3 f^8 - a^6 b c^2 d^3 f^8 - a b^2 c^6 d^3 f^8 -
a^3 b^2 c d^6 f^8 - a b^3 c^2 d^6 f^8 - a^2 b c^3 d^6 f^8)
by substitutions g/. {a^5->-t a^4-p a^3-q a^2-r a-s,b^5->-t b^4-p b^3-q
b^2-r b-s,
c^5->-t c^4-p c^3-q c^2-r c-s,d^5->-t d^4-p d^3-q d^2-r d-s,f^5->-t
f^4-p f^3-q f^2-r f-s}
etc.
How we can do full SymbolicPolynomialMod on the form g after such
procedure will be not higher degree as 4 for a,b,c,d,f ?
Best wishes
Artur
da...@wolfram.com
You will want to use PolynomialReduce in one way or another. You could
simply reduce by what amounts to the rules you indicate above. That does
not really do anything to enforce that the roots are distinct, because
those defining polynomials alone do not enforce distinctness. For that,
one can use the set of relations between roots and coefficients. I'll show
this approach. I use Array to get indexed coefficients and roots, so will
translate your polynomial accordingly.
In[19]:= deg = 5;
coeffs = Array[c, deg, 0];
poly = x^deg + coeffs.x^Range[0, deg - 1];
rts = Array[z, deg];
clist = CoefficientList[Apply[Times, x - rts] - poly, x]
Out[52]= {-c[0] - z[1] z[2] z[3] z[4] z[5], -c[1] +
z[1] z[2] z[3] z[4] + z[1] z[2] z[3] z[5] + z[1] z[2] z[4] z[5] +
z[1] z[3] z[4] z[5] + z[2] z[3] z[4] z[5], -c[2] - z[1] z[2] z[3] -
z[1] z[2] z[4] - z[1] z[3] z[4] - z[2] z[3] z[4] - z[1] z[2] z[5] -
z[1] z[3] z[5] - z[2] z[3] z[5] - z[1] z[4] z[5] - z[2] z[4] z[5] -
z[3] z[4] z[5], -c[3] + z[1] z[2] + z[1] z[3] + z[2] z[3] +
z[1] z[4] + z[2] z[4] + z[3] z[4] + z[1] z[5] + z[2] z[5] +
z[3] z[5] + z[4] z[5], -c[4] - z[1] - z[2] - z[3] - z[4] - z[5]}
That gives the set of defining relations between roots and coefficients.
I'll form a Groebner basis for use in later reduction.
In[25]:= Timing[
gb = GroebnerBasis[clist, rts,
MonomialOrder -> DegreeReverseLexicographic,
CoefficientDomain -> RationalFunctions];]
Out[25]= {0.016, Null}
Fast. Turns out it is the reduction that is slow. Next we translate from
your notation to mine.
g2 = g /. Thread[{a, b, c, d, f} -> rts];
In[30]:= LeafCount[g2]
Out[30]= 2221
We now do the reduction, using the same term ordering as before. Might
want to get a cup of coffee...
In[42]:= Timing[{mults2, red2} =
PolynomialReduce[g2, gb, rts,
MonomialOrder -> DegreeReverseLexicographic,
CoefficientDomain -> RationalFunctions];]
Out[42]= {375.339, Null}
In[43]:= LeafCount[red2]
Out[43]= 22638
That's far bigger than I want to look at, let alone put in email.
Daniel Lichtblau
Wolfram Research
Daniel Lichtblau
> 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?
Reduce can be restricted in that way. Solve cannot, but this will change
in the next release. NSolve cannot be so restricted, at least not in a
useful way. We might try to add such functionality in future but have
not done so as yet.
>
> 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?
No. And I suspect any "global" method will run afoul of complexity, even
when equations are "only" quadratic. But see remark below re optimization.
> 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?
This part I do not understand. The fact that solutions represent triples
in R^3 does not strike me as being relevant to solvers (except insofar
as you impose realness, which might be used to avoid complex space).
> 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).
If you can get by with just one solution (e.g. using symmetries to
generate others), there might be ways to set it up as a constrained
minimization and invoke FindMinimum with interior point or quadratic
programming option. This is not very well documented functionality but
some might be found in the advanced tutorial material for optimization.
> Thanks,
> Robert Hoy
> [...]
Sorry I cannot be of much help here. I simply am not aware of
technologies, other than perhaps quadratic programming, that will find
all solutions (or all in some restricted region) and that might behave
nicely for such problems.
Daniel Lichtblau
Wolfram Research
Robert Hoy
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
DrMajorBob
Ray Koopman
distances instead of the coordinates. If dd is the matrix of squared
distances, and h = IdentityMatrix[n]-1/n, then -.5(h.dd.h) should have
3 positive and n-3 approximate-zero eigenvalues. Sort the eigenvalues,
drop the last 3, minimize the sum of the absolute values of the rest
with respect to the unknown distances. Constraints, such as some
distances being known multiples of others, are easy to impose.
To get the coordinates (with the origin at the centroid of the
configuration), multiply the unit-length eigenvectors corresponding
to the nonzero eigenvalues by the square roots of their eigenvalues.
da...@wolfram.com
I looked again at the original question and my first response. I think it
was not off the mark but I should point out some refinements that might
improve the reach in terms of number of atoms/particles/degrees of
freedom.
The idea is to realize that we can do a bit more to remove degrees of
freedom, by fixing a few coordinates (as you did), and then presolving for
a small subset configuration. For example, since first, second, and fourth
atoms are mutually one unit apart, we know they comprise an equilateral
triangle. We can place it in the x-y coordinate plane with one at origin,
another one unit away on positive x axis, third above x=1/2. This by no
means removes all degrees of freedom, but it provides a nice start.
Here is some code I used.
vars = Array[x, {6, 3}];
amat = {{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}};
subs = Join[Thread[vars[[1]] -> {0, 0, 0}],
Thread[vars[[2]] -> {1, 0, 0}],
Thread[vars[[4]] -> {1/2, Sqrt[3]/2, 0}]];
polys = DeleteCases[
Union[Expand[
Map[Norm[vars[[#[[1]]]] - vars[[#[[2]]]]]^2 - 1 &,
Position[amat, 1]] /. subs /. Abs -> Identity]], 0];
Now the solving step for the remaining polynomials is significantly lower.
In[175]:= Timing[
solns = NSolve[polys + RandomReal[.00001, Length[polys]]];]
Out[175]= {0.296, Null}
Notice I again needed to inject a bit of random perturbation. This is
because I obtained an infinite solution set otherwise (of one dimension).
Offhand I do not know why the system degenerates in this case. That is to
say, I expected that it would be combinatorially rigid, and was clearly
mistaken. In slightly more detail, a GroebnerBasis computation shows that
the sixth atom has x coordinate fixed at 1/2, and the y/z pair lie
anywhere on a circle of radius Sqrt[3]/2 with center on the x axis. So
there is some twisting allowed. Maybe I'll try to model this with
toothpicks, or Mathematica's Manipulate. (Maybe I'll sprout wings and fly
away.)
Anyway, to move onward. The next bit of code removes solutions with
nontrivial imaginary parts.
vars2 = Variables[polys];
s2 = Chop[vars2 /. solns, 10^(-4)];
keep = Map[Apply[And, Map[Head[#] =!= Complex &, #]] &, s2];
realsols = Pick[solns, keep, True];
Now we bring in the inequalities to cull the final viable solutions.
In[187]:=
ineqs = Map[# >= .98 &,
DeleteCases[
Union[Expand[
Map[Norm[vars[[#[[1]]]] - vars[[#[[2]]]]]^2 &,
Position[amat, 0]] /. subs /. Abs -> Identity]], 0]];
In[189]:= goodsols = Pick[realsols, Apply[And, ineqs /. realsols, 1]]
Out[189]= {{x[3, 3] -> -0.816496, x[6, 2] -> -0.288676,
x[5, 3] -> -0.816493, x[6, 3] -> -0.816495,
x[5, 1] -> -1.68335*10^-6, x[6, 1] -> 0.500002, x[3, 2] -> 0.577347,
x[5, 2] -> 0.577347, x[3, 1] -> 0.999997}, {x[3, 3] -> 0.816496,
x[6, 2] -> -0.288676, x[5, 3] -> 0.816493, x[6, 3] -> 0.816495,
x[5, 1] -> -1.68335*10^-6, x[6, 1] -> 0.500002, x[3, 2] -> 0.577347,
x[5, 2] -> 0.577347, x[3, 1] -> 0.999997}}
I think this recovers your octahedron configuration. One solution
configuration lies above the x-y axis, the other below. I do not know of a
convenient way to impose "aboveness" (say) in advance, and still use
methods that avoid computational real algebraic geometry e.g. cylindrical
decomposition, which can be quite expensive.
The upshot is that it is helpful to incorporate early in the process, to
whatever extent possible, removal of degrees of freedom. This tends to
make the hard core part of the solving both faster and more reliable.
If you have not done so already, you might want to check the literature on
methods for solving such systems. These tend to arise in the context of
constraint geometry. There is a fair amount of work from solid modelling
and physical chemistry research (adjust terminology as needed for
correctness). I believe methods involve some graph theory (for
determination of generic combinatorial rigidity), polynomial solving via
various means (including the likes of Mathematica's NSolve), and more. A
place to look for references might be past ADG workshops (Automated
Deduction in Geometry). I recall sure this was a big topic at ADG 2004,
for instance.
Daniel Lichtblau
Wolfram Research