MathGroup Archive 1996

[Date Index] [Thread Index] [Author Index]

Search the Archive

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



  • Prev by Date: Plotting with a simplex
  • Next by Date: Re: InterpolatingFunction::dmwarn: Warning
  • Previous by thread: Solving non-linear equations
  • Next by thread: If[] and % problem