10

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

J. M.'s missing motivation
  • 124,525
  • 11
  • 401
  • 574
Jerry Guern
  • 4,602
  • 18
  • 47
  • 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}. – march Dec 24 '15 at 03:55
  • @march Well then I'm puzzled. I just cut-pasted the code above in a new session and got the exact same wrong results. – Jerry Guern Dec 24 '15 at 03:58
  • I shut down Mathematica, re-started, and ran it again, and I still got that the code works. (V10.0.1 on a Max OS 10.10.5) – march Dec 24 '15 at 04:01
  • 1
    I get the result shown in the post. windows 7, 10.3.1. screen shot: Mathematica graphics and Mathematica graphics – Nasser Dec 24 '15 at 04:08
  • Thanks for the confirmation, @Nasser. Does this mean this is a bug is MMa? – Jerry Guern Dec 24 '15 at 04:22
  • Just to confirm that either version or OS matters, here's my screenshot showing it's working. – march Dec 24 '15 at 04:37
  • I would think it is a bug. We have two different results on two different versions. One of them must be wrong :) – Nasser Dec 24 '15 at 04:48
  • Oh! I am getting correct zeros when I set the WorkingPrecision to 100 as suggested in the Answer below. – Jerry Guern Dec 24 '15 at 04:56
  • But I think the result, from same code, run on 2 different OS's should be the same. Changing WorkingPrecision -> 100 to make the result the same is really a workaround. For me, same code should give same output. Period. Else I would say it is a bug. But I am not an expert on this. – Nasser Dec 24 '15 at 05:02
  • Different OS's has a different compilation of the MMa built-in functions due to features of their relation with CPU, memory etc. And therefore the MachinePrecision definition is a little bit different on Mac and Win machines.. – Rom38 Dec 24 '15 at 05:26
  • Setting the WorkingPrecision to 10 instead of 100 also fixes the problem. Looking more an more like a bug. Is there some official place for reporting bugs? – Jerry Guern Dec 24 '15 at 07:09
  • To send bug report, email support@wolfram.com, make sure you have self contained example of the problem in there and system specs. – Nasser Dec 24 '15 at 09:24
  • 3
    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. – Daniel Lichtblau Dec 24 '15 at 20:01
  • @DanielLichtblau Er… Why it's from a bug if Method->"CompanionMatrix" does better? – xzczd Dec 29 '15 at 03:08
  • 2
    @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). – Daniel Lichtblau Dec 29 '15 at 15:52

3 Answers3

8

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}} *)

Rom38
  • 5,129
  • 13
  • 28
  • This is "only" a 12th order polynomial, and setting the WorkingPrecision to 10 fixes the problem too. I think this is a bug. – Jerry Guern Dec 24 '15 at 07:08
  • 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. – Rom38 Dec 24 '15 at 08:15
  • 2
    @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 – xzczd Dec 25 '15 at 02:14
7

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

John McGee
  • 2,538
  • 11
  • 15
  • 1
    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]. – Jack LaVigne Dec 24 '15 at 21:06
  • 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. – Jerry Guern Dec 25 '15 at 05:26
  • 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 – Jerry Guern Dec 25 '15 at 07:03
  • @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. – Jerry Guern Dec 26 '15 at 02:45
  • @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? – John McGee Dec 26 '15 at 13:35
  • @JohnMcGee I will have to solve this for q1 = PolynomialQuotient[pp[r] - p, pp[1] - p, p]; q2 = PolynomialRemainder[D[pp[r], p] - 1, q1, p]; where r = any prime number. I gave up waiting for r=17, and I want to go as high as possible. BTW, I only need the polynomial in c, I want to eliminate the p's completely. – Jerry Guern Dec 26 '15 at 19:03
4

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}  *)
Michael E2
  • 235,386
  • 17
  • 334
  • 747
  • 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. – Jerry Guern Dec 25 '15 at 05:21
  • @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. – Michael E2 Dec 25 '15 at 05:54