Services & ResourcesWolfram ForumsMathGroup Archive

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

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:
  • Solve problems
    • From: Joerg Schaber <schaber@molgen.mpg.de>