Checking if the roots of a function are real

Question

I'm trying to determine if the roots of a function are real. How would you do that?

(In particular I'm interested in verifying that the roots of LegendreP[6, x] are real.)

Answer 1

This question is not trivial as it would seem and a detailed discussion could help to understand the issue, especially when we deal with roots of special functions, however to do the task as simply as possible this would be the best way :

f[x_] := LegendreP[6, x] 
Reduce[f[x] == 0, x, Reals] == Reduce[f[x] == 0, x]

True

Reduce[f[x] == 0, x, Reals] selects only real roots of f[x] i.e. LegendreP[6, x], while Reduce[f[x] == 0, x] finds all roots. The equality means that we have all roots real.

This way is good enough, because one can easily verify that LegendreP[6, x] is a real function, but if it hadn't been the case we should have rather used another condition assuring the expression to be real, e.g. assume we have another function h[x] not necessarily real and we'd like to find when f[x] == h[x] under assumption that f[x] and x are real, then we could do this :

Reduce[ f[x] == h[x] && Element[ f[x], Reals], x, Reals]

Edit

There are many ways to find only real roots of a polynomial

PolynomialQ[ f[x], x]

True

however a method above based on Reduce is the simplest and powerful enough to verify that we have all roots real, just because it yields the results in the form Root objects. Another way would be Root, but now to write all roots one should know that we have 6 roots (it is not difficult since Exponent[f[x], x] yields 6).

Root[ f[x], #] & /@ Range[6]

{ Root[-5 + 105 #1^2 - 315 #1^4 + 231 #1^6 &, 1],  
  Root[-5 + 105 #1^2 - 315 #1^4 + 231 #1^6 &, 2], 
  Root[-5 + 105 #1^2 - 315 #1^4 + 231 #1^6 &, 3], 
  Root[-5 + 105 #1^2 - 315 #1^4 + 231 #1^6 &, 4], 
  Root[-5 + 105 #1^2 - 315 #1^4 + 231 #1^6 &, 5], 
  Root[-5 + 105 #1^2 - 315 #1^4 + 231 #1^6 &, 6]}

Selecting only real roots yields the same number of roots

Select[%, Re[#] == # &] // Length

6

To elaborate a bit, let's try the other methods e.g. Roots, Solve . To write all root in lists we can do the following :

r = List @@ Roots[f[x] == 0, x][[All, 2]];
s = Solve[f[x] == 0, x][[All, 1, 2]];

here we show e.g. s :

As one can see the output does not seem apparently real just as in case of writing the output of Reduce in the form of radicals :

ToRadicals @ Reduce[ f[x] == 0, x]

One way to assure that roots are only real would be applying FullSimplify with an appropriate ComplexityFunction to radicals, but it can take too much time, e.g. try FullSimplify @ s. Instead of simplifying algebrically the expression, one can check the numerical result, e.g.

Im @ N[s, 100]

{0.*10^-101, 0.*10^-101, 0.*10^-101, 0.*10^-101, 0.*10^-101,  0.*10^-101}

we can apply N directly to Reduce as well :

N @ Reduce[ f[x] == 0, x]

x == -0.238619 || x == 0.238619 || x == -0.661209 || 
x ==  0.661209 || x == -0.93247 || x == 0.93247

If we assume directly only real solutions to be found with Solve, they will be represented by Root objects as in the case of Reduce :

Solve[ f[x] == 0, x, Reals][[All, 1, 2]] == List @@ Reduce[ f[x] == 0, x, Reals][[All, 2]]

True

Now it'll be enlightening to plot the roots :

roots = Tuples[{ List @@ Reduce[ f[x] == 0, x][[All, 2]], {0}}];
Plot[ f[x], {x, -1, 1},
            PlotStyle -> Thick, 
            Epilog -> {PointSize[0.02], Darker @ Green, Point @ roots}]

and to show the relevant structure of real and imaginary parts of LegendreP[6, x] in the complex plane :

GraphicsGrid[{
    { Plot3D[ Re @ f[x + I y], {x, -1, 1}, {y, -1, 1}, Filling -> 0, 
              FillingStyle -> {{Opacity[0.5], Lighter @ Blue}, {Opacity[0.5], Green}}, 
              PlotPoints -> 100, MaxRecursion -> 4 ], 
      Plot3D[ Im @ f[x + I y], {x, -1, 1}, {y, -1, 1}, Filling -> 0, 
              FillingStyle -> {{Opacity[0.5], Lighter@Blue}, {Opacity[0.5], Green}}, 
              PlotPoints -> 100, MaxRecursion -> 4]}, 

    { ContourPlot[ Re @ f[x + I y], {x, -1, 1}, {y, -1, 1}, 
                   Epilog -> {PointSize[0.02], Darker[Green, 0.7], Point @ roots}], 

      ContourPlot[ Im @ f[x + I y], {x, -1, 1}, {y, -1, 1}, 
                   Epilog -> {PointSize[0.02], Darker[Green, 0.7], Point @ roots}]}}]

The option Filling -> 0 helps us to mark the level Re @ f[z] == 0 and Im @ f[z] == 0 in Plot3D's and the green points on ContourPlot's denote the real roots in the complex plane as in the Plot above of the function in the real domain.

Answer 2

I might as well. There is a method that works only for polynomials, but this works for checking if all the roots of a polynomial are real. The method is based on work by Miroslav Fiedler and Gerhard Schmeisser. Briefly, the method constructs (or tries to construct) a tridiagonal companion matrix from your polynomial (i.e. a tridiagonal matrix whose characteristic polynomial is your given polynomial) using a modification of the Euclidean algorithm. On the other hand, even if you do not want the eigenvalues of this special tridiagonal matrix, the entries afford some diagnostic capabilities for checking if your polynomial has all its roots real. Here's how to build the diagonal and one of the off-diagonals:

p = LegendreP[6, x]; n = Exponent[p, x];
p0 = p/Coefficient[p, x, n];
p1 = D[p0/n, x];
{d, e} = MapAt[Most, Transpose[Reap[
     Do[
      {q, r} = PolynomialQuotientRemainder[p0, p1, x];
      s = If[TrueQ[r == 0], 0, Coefficient[-r, x, Exponent[r, x]]];
      If[k < n, p0 = p1;
       p1 = If[TrueQ[r == 0], (#/Coefficient[#, x, Exponent[#, x]]) &[D[p1, x]], -r/s]];
      Sow[{-Coefficient[q, x, 0], s}];
      , {k, n}]][[2, 1]]], 2]
{{0, 0, 0, 0, 0, 0}, {5/11, 8/33, 5/21, 8/35, 1/5}}

The result needed is this: if all the entries of the off-diagonal e are nonnegative, then all the roots of the polynomial are real. (See the papers linked to for more details.) Clearly, this is the case:

And @@ NonNegative[e]
True

If you want to see the actual roots themselves,

Eigenvalues[N[SparseArray[{Band[{1, 1}] -> d,
     Band[{2, 1}] -> Sqrt[e], Band[{1, 2}] -> Sqrt[e]}], 20]] // Sort
{-0.93246951420315202781, -0.66120938646626451366,
-0.23861918608319690863, 0.23861918608319690863,
0.66120938646626451366, 0.93246951420315202781}

or alternatively, Eigenvalues[N[SparseArray[{Band[{1, 1}] -> d, Band[{2, 1}] -> e, Band[{1, 2}] -> 1}], 20]] // Sort.

Compare:

x /. NSolve[LegendreP[6, x], x, 20]
{-0.93246951420315202781, -0.66120938646626451366,
-0.23861918608319690863, 0.23861918608319690863,
0.66120938646626451366, 0.93246951420315202781}

Answer 3

If you are only interested in roots of LegendreP, then you might be better off using a theorem (found e.g. in Szegő, Orthogonal Polynomials) which states:

The zeros of the orthogonal polynomials $p_n(x)$, associated with the distribution $d\alpha(x)$ on the interval $[a, b]$, are real and distinct and are located in the interior of the interval $[a, b]$.

The part of its proof you are interested in - which is the fact that all roots are real - makes use of Sturmian sequences, so perhaps you might want to try playing around with them.