Re: Trouble with coupled quadratic equations where the
- To: mathgroup at smc.vnet.net
- Subject: [mg108366] Re: [mg106476] Trouble with coupled quadratic equations where the
- From: danl at wolfram.com
- Date: Mon, 15 Mar 2010 00:07:05 -0500 (EST)
- References: <201001141046.FAA19729@smc.vnet.net>
> 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
> [...]
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