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 [email protected]

Improve this question
$\endgroup$
15
  • $\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, 2015 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, 2015 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, 2015 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, 2015 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, 2015 at 15:52

3 Answers 3

Reset to default
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}} *)

Improve this answer
$\endgroup$
3
  • $\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, 2015 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, 2015 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, 2015 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

Improve this answer
$\endgroup$
6
  • 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, 2015 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, 2015 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, 2015 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, 2015 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, 2015 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}  *)
Improve this answer
$\endgroup$
2
  • $\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, 2015 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, 2015 at 5:54

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge that you have read and understand our privacy policy and code of conduct.

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