Re: Solving non-linear equations
- To: mathgroup at smc.vnet.net
- Subject: [mg5514] Re: Solving non-linear equations
- From: danl (Daniel Lichtblau)
- Date: Sat, 14 Dec 1996 19:26:10 -0500
- Organization: Wolfram Research, Inc.
- Sender: owner-wri-mathgroup at wolfram.com
In article <58at18$o1v at dragonfly.wolfram.com> elisha at dot.net.au (Luci
Ellis) writes:
> Hello Mathgroup.
>
> One of my colleagues and I are trying to solve out first-order
conditions
> from an optimisation problem. I had written a neat little function that
> does all the necessary transformations. We can get the first order
> conditions really easily, but they seems to eat up so much RAM that it
just
> refuses to solve -- instead the kernel quits.
>
> Here is a simplified version of the problem. Even it doesn't work and we
> REALLY want to do a more complex version.
>
> It's okay if there are multiple solutions, we are expecting that. We
want
> to solve it out analytically to ensure that we have global minima, so
> FindRoot isn't the answer.
>
> Anyway, we want to solve this problem
>
> focs = {0.0175051*a1^2*a2*a3 + 0.25669*a1*a2*a3*d1 - 193.*d2*d3 +
>
> 0.0175051*a1*d1*d2*d3 + 0.25669*d1^2*d2*d3 == 0,
>
> -193.*a2*a3 + 0.0145715*a1^2*a2*a3 + 0.0175051*a1*a2*a3*d1 +
>
> 0.0145715*a1*d1*d2*d3 + 0.0175051*d1^2*d2*d3 == 0,
>
> -0.26055*a1*a2^2*a3 + 0.0145715*a1*a2*a3*d2 - 193.*d1*d3 -
>
> 0.26055*a2*d1*d2*d3 + 0.0145715*d1*d2^2*d3 == 0,
>
> -193.*a1*a3 + 33.6978*a1*a2^2*a3 - 0.26055*a1*a2*a3*d2 +
>
> 33.6978*a2*d1*d2*d3 - 0.26055*d1*d2^2*d3 == 0,
>
> -193.*a1*a2 + 0.25669*a1*a2*a3^2 - 0.42074*a1*a2*a3*d3 +
>
> 0.25669*a3*d1*d2*d3 - 0.42074*d1*d2*d3^2 == 0,
>
> -0.42074*a1*a2*a3^2 - 193.*d1*d2 + 33.6978*a1*a2*a3*d3 -
>
> 0.42074*a3*d1*d2*d3 + 33.6978*d1*d2*d3^2 == 0}
>
> with respect to the algebraic objects: vars = {d1, a1, d2, a2, a3, d3}
>
> And thus get replacement rules {d1 -> some real number, d2 -> another
real
> number... and so on.
>
> Does anyone know why Solve[focs,vars] eats so much RAM?
>
> Please email me at elisha at dot.net.au
>
> Thanks in anticipation.
>
> Luci Ellis
>
>
The code below finds many approximate solutions. I'll indicate some of the
wheres and whyfores along the way.
polys = {0.0175051*a1^2*a2*a3 + 0.25669*a1*a2*a3*d1 - 193.*d2*d3 +
0.0175051*a1*d1*d2*d3 + 0.25669*d1^2*d2*d3,
-193.*a2*a3 + 0.0145715*a1^2*a2*a3 + 0.0175051*a1*a2*a3*d1 +
0.0145715*a1*d1*d2*d3 + 0.0175051*d1^2*d2*d3,
-0.26055*a1*a2^2*a3 + 0.0145715*a1*a2*a3*d2 - 193.*d1*d3 -
0.26055*a2*d1*d2*d3 + 0.0145715*d1*d2^2*d3,
-193.*a1*a3 + 33.6978*a1*a2^2*a3 - 0.26055*a1*a2*a3*d2 +
33.6978*a2*d1*d2*d3 - 0.26055*d1*d2^2*d3,
-193.*a1*a2 + 0.25669*a1*a2*a3^2 - 0.42074*a1*a2*a3*d3 +
0.25669*a3*d1*d2*d3 - 0.42074*d1*d2*d3^2,
-0.42074*a1*a2*a3^2 - 193.*d1*d2 + 33.6978*a1*a2*a3*d3 -
0.42074*a3*d1*d2*d3 + 33.6978*d1*d2*d3^2};
vars = {d1,a1,d2,a2,a3,d3};
First we want to manipulate the polynomials algebraically. It will help to
have exact coefficients. We can use NSolve later to return to approximate
root-finding technology, when we want it.
rpolys = Rationalize[polys, 0];
If we factor them (to break it into smaller subproblems) we find that they
are irreducible. So we go to the next step, and form a Groebner basis
(this will require version 3.0, I think; I doubt previous versions will
complete the task). It turns out that this will now contain polynomials
that factor.
gb = GroebnerBasis[rpolys, vars];
Now we will factor and remove multiplicity; that is, factors of the form
factor^e will be transformed simply to factor. This will not alter the
solution set, though it changes the multiplicity.
expr2 = Map[FactorList, gb];
expr3 = Map[#[[1]]&, expr2, {2}];
expr4 = expr3 /. List[coef_?NumberQ, x___] -> List[x];
If you print out the list, you find that many factors are singletons, that
is, things like a2 or d3. We remove these. This may have the effect of
destroying solutions, but it is necessary because if we try to call Solve
with this processed input as is, it will run out of space attempting to do
LogicalExpand. One can later try to find solutions that involve, say,
a2==0, just by making that substitution and working with the reduced
system. I do not do this below.
LogicalExpand, reform systems of equations, and finally call Solve.
logexpr = Apply[And, Map[Apply[Or,#]&, expr4]];
trunclogexpr = logexpr //. Or[x___,y_Symbol,z___] -> Or[x,z];
logexpr2 = LogicalExpand[trunclogexpr];
logexpr3 = logexpr2 /. And[x___]:>Thread[List[x]==0]
logexpr3 = logexpr2 /. And[x___]:>List[x]
logexpr3 = logexpr2 /. And->List
logexpr4 = Map[Thread[#==0]&, logexpr3]
logexpr5 = logexpr4 /. List->And
nsoln = NSolve[logexpr5, vars]
We now have all results that were not lost from pruning.
We check them all below.
polys /. nsolns // Chop[#, 10^-8]&
I'll now show one case where we explicitly set a variable equal to zero,
and look for solutions.
In[12]:= GroebnerBasis[rpolys /. a2->0, vars] // InputForm
Out[12]//InputForm= {d2*d3, a1*a3, d1*d3, d1*d2}
So if a2 is zero, we have solutions where three other variables are zero
(to see this, just do Solve[%==0, vars]), and the remaining variables are
undetermined. that is, the solution set has positive dimension for these
components.
Daniel Lichtblau
Wolfram Research
danl at wolfram.com