Mathematica 9 is now available
Services & Resources / Wolfram Forums
-----
 /
MathGroup Archive
2006
*January
*February
*March
*April
*May
*June
*July
*August
*September
*October
*November
*December
*Archive Index
*Ask about this page
*Print this page
*Give us feedback
*Sign up for the Wolfram Insider

MathGroup Archive 2006

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

Search the Archive

Re: Solve problems

  • To: mathgroup at smc.vnet.net
  • Subject: [mg64247] Re: [mg64233] Solve problems
  • From: Andrzej Kozlowski <akoz at mimuw.edu.pl>
  • Date: Wed, 8 Feb 2006 03:53:34 -0500 (EST)
  • References: <200602070836.DAA29884@smc.vnet.net>
  • Sender: owner-wri-mathgroup at wolfram.com

On 7 Feb 2006, at 08:36, Joerg Schaber wrote:

> Hi,
>
> I have a system of polynomial equations where Solve cannot find the
> right solutions. Any hints? By the way, actually I just want to  
> find the
> steady states of a differential equation system. If there is another
> clever way, please let me know.
> There are also coonstraints that all variable must be >=0, but  
> including
> those constraints and using Reduce also does not yield a valid  
> solution.
> There exisits a valid solution. I checked this solving the  
> differential
> equation system with NDSolve and let it run into the steady state. But
> this is not very elegant and principally the steady states can be
> calculated directly.
>
> eqns1={0 == 0. c1 + 0.001 c3 - 721.9 c1 c4 + 0.001 c2 c8,
> 0 == 0.001 c3 - 721.9 c1 c4 + 0.001 c5 - 346.09 c4 c6,
> 0 == 0.001 c5 - 346.09 c4 c6 - 989.77 c6 c7 + 0.001 c8,
> 0 == -989.77 c6 c7 + 0.001 c8, c1 + c2 + c3 == 5.7,
> c3 + c4 + c5 == 19.3,
> c5 + c6 + c8 == 4.,
> c7 + c8 == 1.};
>
> sol = Solve[eqns1, {c1, c2, c3, c4, c5, c6, c7, c8},  
> VerifySolutions ->
> True];
>
> \!\({{c2 -> 0.`, c3 ->
>          5.69999921807596`, c5 -> 3.502051416006559`, c1 -> \
> 7.819240400365867`*^-7, c7 -> 0.5020524180816763`, c8 ->
> 0.4979475819183236`, \
> c6 -> 1.0020751177226038`*^-6, c4 -> 10.097949365917481`}, {c2 -> 0.`,
> c3 -> \
> 5.699999640315395`,
>        c5 -> \(-8.352092810365317`\), c1 -> 3.5968460447115093`*^-7,
>           c7 -> \(-11.352093909700997`\), c8 ->  
> 12.352093909700997`, c6 -> \
> \(-1.0993356796497231`*^-6\), c4 -> 21.95209317004992`}, {c2 -> 0.`,
>        c3 -> 8.687577015258206`, c5 -> 10.612427012862872`,
>        c1 -> \(-2.987577015258207`\), c7 -> \ 
> (-1.3272190975834564`*^-7\),
>          c8 -> 1.0000001327219097`, c6 -> \(-7.612427145584781`\),
>         c4 -> \(-4.028121078693263`*^-6\)}, {c2 -> 0.`, c3 -> \
> 19.300001049304143`,
>          c5 -> 9.165049916617935`*^-7, c1 -> \(-13.600001049304144`\),
> c7 -> \
> \(-3.0000004306092727`\), c8 -> 4.000000430609273`,
>          c6 -> \(-1.3471142644128086`*^-6\), c4 -> \
> \(-1.965809135229149`*^-6\)}, {c2 -> 0.`,
>          c3 -> 24.324407481766347`, c5 -> \(-5.024405672582244`\),  
> c1 -> \
> \(-18.624407481766347`\), c7 -> 1.259078407457304`*^-7, c8 -> \
> 0.9999998740921593`, c6 -> 8.024405798490085`,
>         c4 -> \(-1.8091841042296682`*^-6\)}}\)
>
>
> Einsetzen ergibt:
>
> eqns1 /. sol
>
> {{True, False, False, False, True, True, True, True}, {True, True,
>      False, False, True, True, True, True}, {True, False, False, True,
>      True, True, True, True}, {True, False, True, True, True, True,
>      True, True}, {False, False, False, True, True, True, True, True}}
>
>
> best wishes,
>
> joerg
>


This looks to me like a problem with numerical ill conditioning. In  
such cases there is little one can do with MachinePrecision  
arithmetic (at least without subjecting the equations to careful  
numerical analysis). However, one thing one can do is Rationalize and  
use GroebnerBasis. In fact numerical GroebnerBasis with fairly high  
precision usually works fine (and faster than exact GroebnerBasis).  
Here is how to get a solution:

eqns1 = {0 == 0. c1 + 0.001 c3 - 721.9 c1 c4 + 0.001 c2 c8,
0 == 0.001 c3 - 721.9 c1 c4 + 0.001 c5 - 346.09 c4 c6,
0 == 0.001 c5 - 346.09 c4 c6 - 989.77 c6 c7 + 0.001 c8,
0 == -989.77 c6 c7 + 0.001 c8, c1 + c2 + c3 == 5.7,
c3 + c4 + c5 == 19.3,
c5 + c6 + c8 == 4.,
c7 + c8 == 1.};



eqs = Apply[Subtract, eqns1, {1}]


{721.9*c4*c1 + 0.*c1 - 0.001*c3 - 0.001*c2*c8,
   -0.001*c3 + 721.9*c1*c4 - 0.001*c5 + 346.09*c4*c6,
   -0.001*c5 + 346.09*c4*c6 + 989.77*c6*c7 - 0.001*c8,
   989.77*c6*c7 - 0.001*c8, c1 + c2 + c3 - 5.7,
   c3 + c4 + c5 - 19.3, c5 + c6 + c8 - 4., c7 + c8 - 1.}


gr = GroebnerBasis[Rationalize[eqs], {c1, c2, c3, c4, c5, c6, c7, c8},
       CoefficientDomain -> InexactNumbers[300]];


sols = NSolve[gr == 0, {c1, c2, c3, c4, c5, c6, c7, c8}];


N[Select[sols, And @@ (eqns1 /. #1) & ]]


{{c1 -> 3.5968460447115093*^-7, c2 -> 0.,
    c3 -> 5.699999640315395, c4 -> 21.95209317004992,
    c5 -> -8.352092810365317,
    c6 -> -1.0993356796497231*^-6,
    c7 -> -11.352093909700997, c8 -> 12.352093909700997}}

Veryfying again:


eqns1 /. %


{{True, True, True, True, True, True, True, True}}

I am sure Daniel Lichtblau can explain the whole situation much better.

Andrzej Kozlowski


  • References:
  • Prev by Date: Re: Solve problems
  • Next by Date: Re: Using Map with a function of more than 1 argument
  • Previous by thread: Re: Solve problems
  • Next by thread: Re: Solve problems