2
$\begingroup$

Coming from this post: What would be the roots of the derivate of this polynom I did a quick check with Mathematica and arrived at a wrong result. I was initially confused but chose to trust Mathematica and think later.

f[x_, n_] := Product[(x - k), {k, 1, n}]
NSolve[D[f[x, 100], x] == 0, x]
NSolve[D[f[x, 101], x] == 0, x]
NSolve[D[f[x, 102], x] == 0, x]

For $n=100$ I get 99 different results, for $n=102$ I get 101 different results, and for $n=101$ I get 100 times $51.$ as a result.

What is going wrong here?

$\endgroup$
  • 4
    $\begingroup$ using Reals as the domain (NSolve[D[f[x,101], x] == 0, x, Reals]) gets rid of the issue. $\endgroup$ – kglr Jul 16 at 8:56
  • $\begingroup$ btw, this issue does not arise in version 9 (windows 10) $\endgroup$ – kglr Jul 16 at 9:04
  • 3
    $\begingroup$ From a comment in the linked question, setting the Precision to 20 also solves the problem. Still, I find this very curious $\endgroup$ – infinitezero Jul 16 at 11:49
  • 2
    $\begingroup$ It is not the best conditioned of problems. $\endgroup$ – Daniel Lichtblau Jul 16 at 14:24
  • 2
    $\begingroup$ @DanielLichtblau: Indeed. Note that this looks a lot like a derivative of Wilkinson's polynomial, which is a well-known ill-conditioned problem in numerical analysis. $\endgroup$ – Michael Seifert Jul 18 at 13:24
2
$\begingroup$

Actually, solutions for this problem with n greater than approx. 10 are all wrong (no matter if they are different or equal to each other). While for n=7 the max. value of D[f[x, 7], x] on solutions is approx. 10^(-10), the max. value of D[f[x, 13], x] on solutions reaches approx. 3.9, too far from zero. Moreover, a simple analysis shows that, for any given n, there must be (n-1) roots of D[f[x, n], x], that any root belongs to the interval between 1 and n, and that there must be one and only one root in every interval between m and (m+1), m=1,2,...,n-1. However, NSolve[D[f[x, 100], x] == 0, x] gives us three roots between 4 and 5, ten roots between 5 and 6, and sixteen roots greater than 100, e.g. x=282.634 (I believe, numbers can vary depending on version of Mathematica, I use 8.0.4). So, if solutions are wrong, it does not matter whether they are different or equal to each other. Much higher precision computations are required.

$\endgroup$
2
$\begingroup$

The slick way to go about this is to use the logarithmic derivative instead in NSolve[], which has a pretty nice formula:

$$\frac{\mathrm d}{\mathrm dx}\log\left(\prod_{k=1}^n (x-k)\right)=\sum_{k=1}^n\frac1{x-k}$$

Obviously, this has the same roots as the derivative. (This is not true in general, but here the polynomial and its derivative are relatively prime.)

That is to say,

NSolve[Sum[1/(x - k), {k, 1, n}] == 0, x]

Here is a plot of the maximum absolute error for the naive and the log derivative approaches, compared with a solution computed at higher precision:

ListLinePlot[Transpose[Table[Block[{sol, r1, r2},
             sol = Sort[x /. NSolve[D[Product[x - k, {k, 1, n}], x] == 0, x, 
                        WorkingPrecision -> 20]];
             r1 = Sort[x /. NSolve[D[Product[x - k, {k, 1, n}], x] == 0, x]];
             r2 = Sort[x /. NSolve[Sum[1/(x - k), {k, 1, n}] == 0, x]];
             {Norm[sol - r1, ∞], Norm[sol - r2, ∞]}], {n, 3, 50}]],
             DataRange -> {3, 50}, PlotLegends -> {"naive", "log derivative"}, PlotRange -> All]

plot

$\endgroup$
  • 1
    $\begingroup$ On an unrelated note, copying images from Wolfram One is a pain in the bum; I had to make do with taking a screenshot instead. $\endgroup$ – J. M. will be back soon Jul 25 at 13:47

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.