Re: non linear system with 8 equations
- To: mathgroup at smc.vnet.net
- Subject: [mg117263] Re: non linear system with 8 equations
- From: Daniel Lichtblau <danl at wolfram.com>
- Date: Sun, 13 Mar 2011 05:26:31 -0500 (EST)
Ruggero wrote:
> 2011/3/11 Daniel Lichtblau <danl at wolfram.com>:
>> 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