Re: Re: Trouble with coupled quadratic
- To: mathgroup at smc.vnet.net
- Subject: [mg108305] Re: [mg108256] Re: [mg106476] Trouble with coupled quadratic
- From: DrMajorBob <btreat1 at austin.rr.com>
- Date: Sat, 13 Mar 2010 07:56:19 -0500 (EST)
- References: <201001141046.FAA19729@smc.vnet.net>
- Reply-to: drmajorbob at yahoo.com
If you have 3n variables, I don't see why you'd think you're in 3-space.
Bobby
On Fri, 12 Mar 2010 06:07:48 -0600, Robert Hoy <robert.hoy at yale.edu> wrote:
> 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
>
>
--
DrMajorBob at yahoo.com