# Solving set of non-linear simultaneous equations gives PolynomialGCD::lgrexp error - what next?

I'm trying to solve a set of six simultaneous non-linear equations to get a numerical answer using NSolve:

NSolve[{-0.5547005383792517 +
0.8082903768654761/E^(0.08452994616207483*(0.10136907301545368 +
beta1/0.3^alpha)) ==
a, -0.3147005383792517 + 0.669726312259966/

E^(0.08452994616207483*(0.10136907301545368 +
beta1/0.45^alpha)) == a,
-0.13470053837925167 +
0.5658032638058333/
E^(0.08452994616207483*(0.10136907301545368 +
beta1/0.55^alpha)) ==
a, -0.540774588761456 +
0.805711228582561/
E^(0.0779225411238544*(0.10136907301545368 + beta2/0.3^alpha)) ==
a, -0.420774588761456 +
0.7324647532568737/
E^(0.0779225411238544*(0.10136907301545368 + beta2/0.45^alpha)) ==
a, -0.160774588761456 +
0.5737640567178843/
E^(0.0779225411238544*(0.10136907301545368 + beta2/0.55^alpha)) ==
a}, {alpha, beta1, beta2, a}, Reals]


(I'm not very experienced with Mathematica - so I hope that is the correct way to paste code in here so that others can copy into Mathematica and try running it)

I get the following errors:

PolynomialGCD::lrgexp: Exponent is out of bounds for function PolynomialGCD.


This seems to suggest that Mathematica can't deal with the large values used as the exponents in the equations given - but as far as I can see the values aren't particularly large. I do know roughly what sort of ranges the results for each of the variables should have (eg. a should be between 0 and 1) but I can't seem to find a way to specify this extra knowledge as an option to NSolve.

I'm not sure where to progress from here, now that I've got these errors. Does anyone have any ideas? I am taking the equations from a paper that I am trying to reimplement, but I'm pretty sure the equations are how they're meant to be. I may be doing something silly in Mathematica though.

You seem to be trying to solve an overdetermined system? Is it supposed to be alpha1 and alpha2 for the top and bottom three equations? If not, you can get an approximate solution by way of least-squares.

eqns = {
-0.5547 + 0.8083 E^-(0.0845*(0.1014 + beta1/0.30^alpha)) - a,
-0.3147 + 0.6697 E^-(0.0845*(0.1014 + beta1/0.45^alpha)) - a,
-0.1347 + 0.5658 E^-(0.0845*(0.1014 + beta1/0.55^alpha)) - a,
-0.5408 + 0.8057 E^-(0.0779*(0.1014 + beta2/0.30^alpha)) - a,
-0.4208 + 0.7324 E^-(0.0779*(0.1014 + beta2/0.45^alpha)) - a,
-0.1608 + 0.5738 E^-(0.0779*(0.1014 + beta2/0.55^alpha)) - a
};

sol = FindMinimum[Total[eqns ^2], {alpha, a, beta1, beta2}, Method -> "LevenbergMarquardt"]
>> {0.000126838, {alpha -> 0.00366703, a -> 0.810457, beta1 -> -6.23349, beta2 -> -6.7527}}


And to test that solution

eqns /. Last@sol
>> {-0.00490439, 0.000966248, 0.00588526, 0.00458171, 0.000264733, -0.00679357}


It doesn't seem to be a very good fit. I'm mostly concerned your equations may be off.

• Hmmm - thanks. I'll have a deeper look into my equations and try and understand whether they're correct or not. The answers given by the least-squares aren't sensible answers for the scientific application that I'm using this for - which suggests that my equations are wrong. – robintw Feb 5 '13 at 21:03

[Not an answer, just too big to fit in the margin.]

The large exponents come from the solver code rationalizing the powers of E and then invoking Together. In this case it is specifically doing\

Together[6040676245532097866985581585354 -
20344251995931118705100736513935*E^(-917705716247/107099481787092 -
(7444009*2^(-2 + alpha)*3^(-1 - alpha)*5^alpha*beta1)/7338631) +
16856665939485786567230379212400*E^(-917705716247/107099481787092 -
(7444009*3^(-1 - 2*alpha)*4^(-1 + alpha)*5^alpha*beta1)/7338631)]


One can learn this by first doing:

Unprotect[Together];
Together[xxx__] := Null /; (Print[InputForm[together[xxx]]]; False)
`