MathGroup Archive 2012

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

Search the Archive

Re: Can I solve this system of nonlinear equations?

  • To: mathgroup at smc.vnet.net
  • Subject: [mg125267] Re: Can I solve this system of nonlinear equations?
  • From: danl at wolfram.com
  • Date: Sat, 3 Mar 2012 06:53:36 -0500 (EST)
  • Delivered-to: l-mathgroup@mail-archive0.wolfram.com
  • References: <jil5ld$rrm$1@smc.vnet.net>

On Wednesday, February 29, 2012 6:28:29 AM UTC-6, Andy wrote:
> I'm dealing with systems of nonlinear equations that have 8 equations
> and 8 unknowns.  Here's an example:
>
> Solve[{(((c - a)/0.002) - (0.995018769272803 + h*b)) == 0,
>   (((d - b)/0.002) - (0.990074756047929 + h*c)) == 0,
>   (((e - c)/0.002) - (0.985167483257382 + h*d)) == 0,
>   (((f - d)/0.002) - (0.980296479563062 + h*e)) == 0,
>   (((g - e)/0.002) - (0.975461279165159 + h*f)) == 0,
>   (((-1*e + 8*d - 8*b + a)/(12*0.001)) - (0.990074756047929 + h*c)) ==
> 0,
>   (((-1*f + 8*e - 8*c + b)/(12*0.001)) - (0.985167483257382 + h*d)) ==
> 0,
>   (((-1*g + 8*f - 8*d + c)/(12*0.001)) - (0.980296479563062 + h*e)) ==
> 0}, {a, b, c, d, e, f, g, h}]
>
> Whenever I try this, Mathematica 7 just returns the empty set {}.  How
> can I tell if this is unsolvable?  Shouldn't I at least be able to get
> a numerical approximation with NSolve?  I've tried using stochastic
> optimization to get approximate answers but every method gives poor
> results, and that's why I would like to at least approximately solve
> this if possible.  Thanks very much for any help~

It is not solvable. Or at least, a rationalized version has no solution. Here is one way to show this.

polys = Apply[
   Subtract, {(((c - a)/0.002) - (0.995018769272803 + h*b)) ==
     0, (((d - b)/0.002) - (0.990074756047929 + h*c)) ==
     0, (((e - c)/0.002) - (0.985167483257382 + h*d)) ==
     0, (((f - d)/0.002) - (0.980296479563062 + h*e)) ==
     0, (((g - e)/0.002) - (0.975461279165159 + h*f)) ==
     0, (((-1*e + 8*d - 8*b + a)/(12*0.001)) - (0.990074756047929 +
         h*c)) ==
     0, (((-1*f + 8*e - 8*c + b)/(12*0.001)) - (0.985167483257382 +
         h*d)) ==
     0, (((-1*g + 8*f - 8*d + c)/(12*0.001)) - (0.980296479563062 +
         h*e)) == 0}, 1];

p2 = Rationalize[polys, 0];
p3 = Expand[p2];
vars = Variables[p2];

For showing non-solvability, it suffices (by Hilbert's Nullstellensatz, apologies for pulling that bit of math jargon out of a hat) to find a way to write a constant as an algebraic sum of these polynomials. That is to say, we show the ideal defined by those polynomials contains 1. One way to do so uses a bit of code that I posted to MathGroup around a year ago.

http://forums.wolfram.com/mathgroup/archive/2011/Mar/msg00362.html

I'll just copy it here for completeness.

moduleGroebnerBasis[polys_, vars_, cvars_, opts___] :=
 Module[{newpols, rels, len = Length[cvars], gb, j, k, rul},
  rels = Flatten[
    Table[cvars[[j]]*cvars[[k]], {j, len}, {k, j, len}]];
  newpols = Join[polys, rels];
  gb = GroebnerBasis[newpols, Join[cvars, vars], opts];
  rul = Map[(# :> {}) &, rels];
  gb = Flatten[gb /. rul];
  Collect[gb, cvars]]


conversionMatrix[polys_, vars_] :=
 Module[{aa, coords, pmat, len = Length[polys], newpolys, mgb, gb,
   convmat, fvar, rvars}, coords = Array[aa, len + 1];
  fvar = First[coords];
  rvars = Rest[coords];
  pmat = Transpose[Join[{polys}, IdentityMatrix[len]]];
  newpolys = pmat.coords;
  mgb = moduleGroebnerBasis[newpolys, vars, coords];
  gb = mgb /. Join[{fvar -> 1}, Thread[rvars -> 0]] /. 0 :> Sequence[];
  convmat = Select/. fvar -> 0;
  {gb, convmat /.
    Thread[rvars -> Table[UnitVector[len, j], {j, len}]]}]

On your example we proceed as follows.

In[40]:= {gb, cmat} = conversionMatrix[p3, vars]

Out[40]= {{2321109176622379997592951205}, \
{{-2476806404274022328337996647266500,
   19814451234192178626703973178132000,
   29721676851288267940055959767198 h, \
-19814451234192178626703973178132000,
   2476806404274022328337996647266500, \
-14860838425644133970027979883599000, \
-29721676851288267940055959767198 h,
   14860838425644133970027979883599000}}}

We check the result. Specifically we show that the conversion matrix, multiplied by the polynomial system input regarded as a vector, indeed gives a constant.

In[42]:= Expand[cmat.p3]

Out[42]= {2321109176622379997592951205}

So there you have it: no solution for this system.

Daniel Lichtblau
Wolfram Research




  • Prev by Date: Re: Can I solve this system of nonlinear equations?
  • Next by Date: Re: Can I solve this system of nonlinear equations?
  • Previous by thread: Re: Can I solve this system of nonlinear equations?
  • Next by thread: Re: Can I solve this system of nonlinear equations?