I am using NSolve to find roots to a system of polynomials that describe some chemical reactions. However, I am seeing some odd behavior when I restrict the domain.
First, the basic setup:
c1 = 50.0;
c2 = 1.5;
c3 = 0.4;
c4 = 0.5;
c5 = 8.0;
c6 = 3.2;
k1 = 3.7*10^-18;
k2 = 8.3*10^-9;
k3 = 8.7*10^-13;
k4 = 3.2*10^-7;
k5 = 8.6*10^-3;
f1 = 1/2*(c1-c2+c3+2*c4-c5+c6-x1-2*x2-x3+x4+(x3*x4)/(k5*x1));
f2 = 1/(2*k5*x1)*(3*k5*c6*x1-k5*c1*x1+2*k5*c4*x1+k5*c3*x1-k5*c2*x1-k5*c5*x1+k5*x1^2-2*k5*x1*x2-k5*x1*x3+k5*x1*x4-x3*x4);
eqs={
f1*f2 - k1*x1^2 == 0,
f1*x3 - k2*x2*x1 == 0,
f1*(c6 - x3 - x2 - (x3*x4)/(k5*x1)) - k3*x1*x3 == 0,
f2*(c5 - x4 - (x3*x4)/(k5*x1)) - k4*x1*x4 == 0
};
NSolve[eqs,{x1, x2, x3, x4}]
This returns 32 complex roots. However, I know at the onset that there is at least one real-valued root, and exactly one real root with all four variables greater than zero. So, I try it again and use the specification suggested in the documentation to search for Reals only.
NSolve[eqs, {x1, x2, x3, x4}, Reals]
This, however, returns {} for no solutions. However, while playing around with the options I accidentally used "Real" instead of "Reals":
NSolve[eqs, {x1, x2, x3, x4}, Real]
Now, I receive a warning:
NSolve::bdomv: Warning: Real is not a valid domain specification. Assuming it is a variable to eliminate. >>
But I also get 8 real roots! Interestingly, these are the same roots I find if I run the original calculation at any WorkingPrecision that is not MachinePrecision...
NSolve[eqs, {x1, x2, x3, x4}, WorkingPrecision->MachinePrecision]
NSolve[eqs, {x1, x2, x3, x4}, WorkingPrecision->2]
NSolve[eqs, {x1, x2, x3, x4}, WorkingPrecision->20]
Lastly, I also tried restricting the domain to positive reals. This fails to return any result. If I combine this with a finite WorkingPrecision I get the correct root. If I instead specify "Reals" I get nothing, and I use "Real" I get warnings about the domain as well as infinite solutions, but the same root as above when I used Working Precision.
NSolve[Catenate[{eqs, {x1 > 0, x2 > 0, x3 > 0, x4 > 0}}], {x1, x2, x3, x4}]
NSolve[Catenate[{eqs, {x1 > 0, x2 > 0, x3 > 0, x4 > 0}}], {x1, x2, x3, x4}, WorkingPrecision -> 8]
NSolve[Catenate[{eqs, {x1 > 0, x2 > 0, x3 > 0, x4 > 0}}], {x1, x2, x3, x4}, Reals]
NSolve[Catenate[{eqs, {x1 > 0, x2 > 0, x3 > 0, x4 > 0}}], {x1, x2, x3, x4}, Real]
{}
{{x1 -> 49.4605, x2 -> 0.00100605, x3 -> 0.511294, x4 -> 2.21332}, {x1 -> 53.0999, x2 -> 4.30992*10^-13, x3 -> 0.000114669, x4 -> 7.99797}}
{}
NSolve::bdomv: Warning: Real is not a valid domain specification. Assuming it is a variable to eliminate. >>
NSolve::infsolns: Infinite solution set has dimension at least 1. Returning intersection of solutions with -((41688 x1)/65167)-(153968 x2)/195501+(153196 x3)/195501+(185938 x4)/195501+(38650 System`NSolveDump`Y$31197[1])/65167 == 1. >>
{{x1 -> 49.4605, x2 -> 0.00100591, x3 -> 0.511294, x4 -> 2.21332}}
Anyone have insight into these behaviors?
