MathGroup Archive 2011

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

Search the Archive

Re: non linear system with 8 equations

  • To: mathgroup at
  • Subject: [mg117263] Re: non linear system with 8 equations
  • From: Daniel Lichtblau <danl at>
  • Date: Sun, 13 Mar 2011 05:26:31 -0500 (EST)

Ruggero wrote:
> 2011/3/11 Daniel Lichtblau <danl at>:
>> wiso wrote:
>>> I've a non linear system with 8 equation:
>>> solution = Solve[{n == ng + nn,
>>>     na == eag ng + ean nn,
>>>    nb == ebg ng + ebn nn,
>>>    nc == ecg ng + ecn nn,
>>>    nab == kabg eag ebg ng + kabn ean ebn nn,
>>>    nac == kacg eag ecg ng + kacn ean ecn nn,
>>>    nbc == kbcg ebg ecg ng + kbcn ebn ecn nn,
>>>    nabc == kabcg eag ebg ecg ng + kabcn ean ebn ecn nn }, {ng, nn,
>>>    eag, ebg, ecg, ean, ebn, ecn}];
>>> my pc is cogitating since more than one hour. Is there some method to
>>> speed it up? All variables are real, n* are positive integer, e* are
>>> real 0<e*<1. k* are known parameters
>> I doubt Solve will have any chance with this. It would need to compute some
>> flavor of GroebnerBasis, and that seems to be a difficult undertaking for
>> this system.
>> If you provide explicit numeric values for those parameters then NSolve will
>> handle it. Below is some code for this. It takes the polynomials, variables
>> of interest, parameter substitutions, and constraints entered as a list.
>> In[160]:=
>> constrainedSolutions[polys_, vars_, paramsubs_, constraints_] :=
>>  Module[
>>  {solns},
>>  solns = NSolve[polys /. paramsubs, vars];
>>  Pick[solns, Quiet[And @@@ Release[constraints /. solns]]
>>   ]]
>> Here is your example. I changed the constraints so as to find parameter
>> substitutions that might yield a few valid solutions. I show one such run
>> below.
>> In[161]:=
>> eqns = {n == ng + nn, na == eag ng + ean nn, nb == ebg ng + ebn nn,
>>   nc == ecg ng + ecn nn, nab == kabg eag ebg ng + kabn ean ebn nn,
>>   nac == kacg eag ecg ng + kacn ean ecn nn,
>>   nbc == kbcg ebg ecg ng + kbcn ebn ecn nn,
>>   nabc == kabcg eag ebg ecg ng + kabcn ean ebn ecn nn};
>> polys = Apply[Subtract, eqns, 1];
>> vars = {ng, nn, eag, ebg, ecg, ean, ebn, ecn};
>> params = Complement[Variables[polys], vars];
>> paramsubs = Thread[params -> RandomReal[{-10, 10}, Length[params]]];
>> SeedRandom[11111];
>> constraints =
>>  Flatten[{Map[Hold[Head[#] === Real] && # >= 0 &, {ng, nn}],
>>    Map[Hold[Head[#] === Real] && -100 <= # <= 100 &, {eag, ebg, ecg,
>>      ean, ebn, ecn}]}];
>> In[168]:= constrainedSolutions[polys, vars, paramsubs, constraints]
>> Out[168]= {{ng -> 7.30836, nn -> 0.254874, eag -> 0.186697,
>>  ebg -> -0.90457, ecg -> -0.0648619, ean -> 2.44744, ebn -> 0.726247,
>>   ecn -> -2.57797}, {ng -> 4.8318, nn -> 2.73144, eag -> -0.428736,
>>  ebg -> -1.02281, ecg -> -0.46866, ean -> 1.48632, ebn -> -0.543236,
>>  ecn -> 0.414936}, {ng -> 1.00957, nn -> 6.55367, eag -> -0.174131,
>>  ebg -> -3.37666, ecg -> -4.50395, ean -> 0.330202, ebn -> -0.460332,
>>   ecn -> 0.521227}}
>> If you know a priori that your system will have exactly one solution
>> satisfying the constraints, then you might do better using FindRoot with
>> appropriate initial values for the variables. But for this particular system
>> of polynomials NSolve will be fast (again, provided you give numerical
>> values up front for the parameters).
>> Daniel Lichtblau
>> Wolfram Research
> Thank you for the reply (why didn't you post it on the ng?). Well, I
> know I can find a numberical solution very quickly, but I need to find
> a closed formula, even expressed as root of a n degree polynomial. My
> big problem is that every know number is known with an error, as n +/-
> delta n, and I want to propagate every error in the solution.

I do not see a viable way to do that in full. Even if you managed to get 
a symbolic "solution", I fail to see how it would help you to gauge 
error bounds. It would be laden with parametrized Root[] objects. I do 
not see how they might be coerced to readily give up their secrets in 
regards to the ranges tehy might occupy. What I'm trying to say is that 
a closed form is not only most likely intractable to find, but also 
unuseful once you have it.

If your coefficients are approximately known, say to 3-4 places (even 2 
might suffice), then you could assess the conditioning of the solution 
set (that is, how much the roots might vary given the coefficients can 
be anywhere in their tolerance region). One way to do that, probably the 
most reliable and certainly the easiest, would be via statistical 
sampling of the coefficient space. You simply take high precision values 
selected according to whatever distribution you find appropriate. 
Gaussian is common for this, where you center the Gaussian at the center 
value of your tolerance range, and set the variance so that the upper 
end is, say, three standard deviations.

If you start from one solution e.g. the one given by taking high 
precision coefficients to be the centers of their allowable ranges, you 
can usually get away with an extremely simple type of homotopy approach. 
Specifically, call FindRoot for the altered systems, with starting 
points set to be the solutions to your first system.

This should give a quite reasonable idea of what are the ranges in which 
the roots might lie, when your coefficients are only imprecisely known 
but at least have modest precision.

If you do not know your coefficients to a few places you might still try 
a sampling approach like this. It will give an idea of your root ranges, 
but it is more difficult to justify that this idea in any way "fills 
out" the range of possibilities.

Daniel Lichtblau
Wolfram Research

  • Prev by Date: Vectorized molecular dynamics with Mathematica?
  • Next by Date: Re: what's new in 8.0.1?
  • Previous by thread: Re: non linear system with 8 equations
  • Next by thread: Re: non linear system with 8 equations