10
$\begingroup$

OP UPDATE: I received an email from WR on 1-18-2016: "...It does appear that the NSolve function is not behaving properly in this case and I have forwarded an incident report to our developers with the information you provided..." So, WorkingPrecision->Automatic is indeed failing, but the workaround is to manually set any value, then NSolve works fine.

OP UPDATE: Email from WR on 11-21-2016: "In December 2015 you reported an issue with Mathematica wherein NSolve returns bad results. We believe that the issue has been resolved in the current release of Mathematica."

The following code creates two polynomials $q_1$ and $q_2$ in variables $c$ and $p$, then uses NSolve to find roots. But the polynomials don't evaluate to zero at those roots.

pp[n_] := If[n > 0, pp[n - 1]^2 + c, p];
q1 = PolynomialQuotient[pp[4] - p, pp[2] - p, p];
q2 = PolynomialRemainder[D[pp[4], p] - 1, q1, p];
soln = NSolve[{q1 == 0 , q2 == 0}, {c, p}];
Mean[Abs[{q1, q2} /. soln]]

$\{1.07185\times 10^{12},1.12062\times 10^{13}\}$

Is this a bug, or am I doing something wrong? I've tried using Eliminate, Solve, and && between my conditions, but I got either no results or wrong results.

OP EDIT: In testing solutions offered below, I found that setting WorkingPrecision to any value, even 10, makes the results come back near the desired zeros. But if I set WorkingPrecision->Automatic then the same wrong results come back. Other users have gotten the same results I did.

I have reported this issue to support@wolfram.com

$\endgroup$
  • $\begingroup$ When I run your code and evaluate q1 /. soln // Chop and q2 /. soln // Chop, I get from both a list of zeros. From your code exactly, I get {4.17433*10^-13,4.75146*10^-12}. $\endgroup$ – march Dec 24 '15 at 3:55
  • $\begingroup$ @march Well then I'm puzzled. I just cut-pasted the code above in a new session and got the exact same wrong results. $\endgroup$ – Jerry Guern Dec 24 '15 at 3:58
  • 1
    $\begingroup$ I get the result shown in the post. windows 7, 10.3.1. screen shot: !Mathematica graphics and !Mathematica graphics $\endgroup$ – Nasser Dec 24 '15 at 4:08
  • 3
    $\begingroup$ Try 'NSolve[...,Method->"CompanionMatrix"]' and see if that does better. That can also then indicate if this is from a bug or just a precision issue. $\endgroup$ – Daniel Lichtblau Dec 24 '15 at 20:01
  • 2
    $\begingroup$ @xzczd It helps to distinguish between whether it is a problem that requires high precision vs. one that is just not being handled well. This one falls into the latter category (I filed a report to that effect). $\endgroup$ – Daniel Lichtblau Dec 29 '15 at 15:52
8
$\begingroup$

These results are strongly dependent on WorkingPrecision settings for NSolve. The default settings are not enough for polynomials with so high order.

pp[n_] := If[n > 0, pp[n - 1]^2 + c, p];
q1 = PolynomialQuotient[pp[4] - p, pp[2] - p, p];
q2 = PolynomialRemainder[D[pp[4], p] - 1, q1, p];
soln = NSolve[{q1 == 0, q2 == 0}, {c, p}, WorkingPrecision -> 100];
rrr =Chop[{q1, q2} /. soln]

The answer is: (*{{0, 0}, {0, 0}, {0, 0}, {0, 0}, {0, 0}, {0, 0}, {0, 0}, {0, 0}, {0, 0}, {0, 0}, {0, 0}, {0, 0}, {0, 0}, {0, 0}, {0, 0}, {0, 0}, {0, 0}, {0, 0}, {0, 0}, {0, 0}, {0, 0}, {0, 0}, {0, 0}, {0, 0}} *)

$\endgroup$
  • $\begingroup$ This is "only" a 12th order polynomial, and setting the WorkingPrecision to 10 fixes the problem too. I think this is a bug. $\endgroup$ – Jerry Guern Dec 24 '15 at 7:08
  • $\begingroup$ Hm, may be you are right about definition of this as bug, because if I use WorkingPrecision->Automatic then I have non-zero answer. But usage of 14 or 16 instead of Automatic (which is $MachinePrecision==15.9546 in my system) gave me desired zero. $\endgroup$ – Rom38 Dec 24 '15 at 8:15
  • 2
    $\begingroup$ @JerryGuern I'd like to mention that Automatic isn't $MachinePrecision but MachinePrecision. Though MachinePrecision == $MachinePrecision returns True, when they're used as the option value of WorkingPrecision, they're totally different, MachinePrecision will lead to… Er… machine precision calculation, while $MachinePrecision will lead to arbitrary precision calculation. (Not sure if I've used proper terminology.) This issue has been mentioned sporadically in many posts, an example is this $\endgroup$ – xzczd Dec 25 '15 at 2:14
7
$\begingroup$

Here is my approach using Groebner bases.

gb = Collect[GroebnerBasis[{q1, q2}, {p, c}], p]

csol = Solve[First@gb == 0, c]  (* the first basis element is a polynomial in c only *)

peqs = Last@gb /. csol  (* substitute the c solutions *)

psol=Solve[# == 0, p] & /@ peqs

The solutions given in csol and psol are exact, as you can verify

$\endgroup$
  • 1
    $\begingroup$ I think this is a great answer but suggest that you move it to the original question. I also found that using t in place of 1 in q2 provided solutions for c[t]. $\endgroup$ – Jack LaVigne Dec 24 '15 at 21:06
  • $\begingroup$ The original question Jack is referring to is THIS. Your answer fits that question better, since this questions is about a bug I found while testing solutions to that question. $\endgroup$ – Jerry Guern Dec 25 '15 at 5:26
  • $\begingroup$ Aaaand now I had to post another question on the Math SE, because I've never studied Ring Theory and couldn't understand what I found online about Groebner Basis. Here it is $\endgroup$ – Jerry Guern Dec 25 '15 at 7:03
  • $\begingroup$ @johnMcGee, if all I want is csol and want to just eliminate the p's, this the fastest way to do that on MMa this: GroebnerBasis[{q1, q2}, {p, c},p]? It gives me what I want, but it takes a veeeeery long time if the polynomials are over degree 10 or so, and I need to apply this to as high degrees as possible. $\endgroup$ – Jerry Guern Dec 26 '15 at 2:45
  • $\begingroup$ @Jerry Guern gb = GroebnerBasis[{q1, q2}, {p, c}, MonomialOrder -> DegreeReverseLexicographic] runs very quickly on my machine with your given polynomials. Can you post the polynomials that run very slowly? $\endgroup$ – John McGee Dec 26 '15 at 13:35
4
$\begingroup$

If you use any of the other monomial orderings than Lexicographic, the solutions are accurate:

Grid@Table[{mo, 
   soln = NSolve[{q1 == 0, q2 == 0}, {c, p}, Method -> {"MonomialOrder" -> mo}]; 
   Mean[Abs[{q1, q2} /. soln]]},
 {mo, {Lexicographic, DegreeLexicographic, DegreeReverseLexicographic, EliminationOrder}}]

Mathematica graphics

I thought of this because switching c and p in the OP's code changed the result:

soln = NSolve[{q1 == 0, q2 == 0}, {p, c}]; Mean[Abs[{q1, q2} /. soln]]
(*  {1.0912*10^14, 1.17552*10^15}  *)
$\endgroup$
  • $\begingroup$ How do you interpret what you've demonstrated here? BTW, I sent the link to this SE page in my report to support@wolfram,com so hopefully they'll see this result too. $\endgroup$ – Jerry Guern Dec 25 '15 at 5:21
  • $\begingroup$ @JerryGuern My guess is that it's a precision thing that depends on the way the polynomial system is broken down. I think it's going to be hard to track it down exactly (maybe very hard, or maybe not), but the WRI folks will probably have tools and insights that we others don't have. $\endgroup$ – Michael E2 Dec 25 '15 at 5:54

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.