Zeros of high degree polynomials

Question

I am working with Hermite polynomials in Mathematica with the built-in function HermiteH. I want to compute the zeros of the polynomial HermiteH[N,x] for N as large as I can, as I am computing a Gaussian quadrature rule for integration of a smooth function.

I observe that, as N increases, the accuracy of NSolve drops dramatically. For instance,

HermiteH[18, x /. NSolve[HermiteH[18, x] == 0, x, Reals]]

{1., -0.28125, -0.0175781, -0.000732422, 0.000366211, -0.0000686646, 0.0000267029,
 9.53674*10^-7, 2.6226*10^-6, 2.6226*10^-6, 9.53674*10^-7, 0.0000267029, -0.0000686646,
 0.000366211, -0.000732422, -0.0175781, -0.28125, 1.}

You see, Hermite[18, x] evaluates at the first root as 1. Is there any way to improve the accuracy? Is there any built-in function to compute the roots of Hermite polynomials? Thank you for your help.

Why not use Solve instead of NSolve? Or add a WorkingPrecision option? — Carl Woll, Jun 28 '19 at 19:56
Don't use single argument N, as this will use machine numbers, and you will lose lots of precision due to subtractive cancellation. — Carl Woll, Jun 28 '19 at 20:05
Except for the 5th and 14th, roots NSolve seems to have found in each case the number x such that Abs[Hermite[18, x]] is as small as possible in a small neighborhood of x. Keep in mind that floating-point numbers are discrete. When the derivative of a function is large, as it is especially at the extreme roots, the difference between function values at adjacent floating-point numbers can seem large. The error in the residual should be bounded by D[HermiteH[18, x], x] x $MachineEpsilon/2 /. {x -> x0}, where x0 is a root found by NSolve. — Michael E2, Jun 29 '19 at 04:48
Possible duplicate: https://mathematica.stackexchange.com/q/51098/4999 — Michael E2, Jun 29 '19 at 04:50

score 9 · Accepted Answer · answered Jun 28 '19 at 20:04

Don't use machine numbers, as subtractive cancellation will cause enormous precision loss, as is common with high order polynomials.

You can either work with exact results using Solve:

HermiteH[18, x /. Solve[HermiteH[18,x]==0,x,Reals]] //Simplify

{0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0}

Or you can use the WorkingPrecision option:

HermiteH[18, x /. NSolve[HermiteH[18,x]==0, x, Reals, WorkingPrecision->50]]

{0.*10^-33, 0.*10^-35, 0.*10^-36, 0.*10^-37, 0.*10^-38, 0.*10^-38, 0.*10^-39, 0.*10^-39, 0.*10^-40, 0.*10^-40, 0.*10^-39, 0.*10^-39, 0.*10^-38, 0.*10^-38, 0.*10^-37, 0.*10^-36, 0.*10^-35, 0.*10^-33}

You can see that subtractive cancellation causes a precision loss of about 17 digits for the first and last roots.

A quibble: In my view, "subtractive cancellation" typical of polynomials is not the culprit here. One cannot avoid loss of (relative) precision in evaluating $f(x)$ if $f(x^)=0$ and $x$ is close (enough) to $x^$; it does not follow that there is a loss of accuracy. Indeed, on numeric inputs HermiteH[integer, real] uses a stable algorithm that avoids the subtractive precision loss that occurs in evaluating a polynomial in the power basis. Observe the extra loss of accuracy when NSolve is moved outside: HermiteH[18, x] /. NSolve[HermiteH[18, x] == 0, x, Reals, WorkingPrecision -> 50]. — Michael E2, Jun 29 '19 at 16:08

Roman · Answer 2 · 2019-06-29T10:19:53.427

For polynomials like HermiteH, these roots are represented in Mathematica as infinite-precision Root objects. The $k$th root of the $n$th Hermite polynomial is, with infinite precision, represented by

R[n_, k_] := Root[HermiteH[n, #] &, k]

What Carl's use of Solve does is simply to make a list of such Root objects. You can work with these objects analytically (using RootReduce etc.), or you can convert them to numerical.

For example, the 7th root of $H_{18}$ would be

r = R[18, 7]
(*    a root around -1.30...    *)

Numerically:

N[r]
(*    -1.30092    *)
N[r, 100]
(*    -1.300920858389617365666265554392610580218134639661226522772309775882782630084141194539623631652544514    *)

analytic transformations:

r^2 // RootReduce
(*    a root around 1.69...    *)
HermiteH[18, r] // RootReduce
(*    0    *)

score 2 · Answer 3 · answered Nov 19 '19 at 01:36

At least for this case, one can also consider getting the eigenvalues of the Jacobi matrix associated with the Hermite polynomials. Recall that these matrices are constructed from the coefficients of the three-term recurrences that generate the corresponding orthogonal polynomial. Applied to this case, we have:

With[{n = 18}, 
     s1 = Sort[Eigenvalues[N[SparseArray[{{j_, k_} /; Abs[j - k] == 1 :> Sqrt[Min[j, k]/2]},
                                         {n, n}]]]]]
   {-5.04836, -4.24812, -3.57377, -2.96138, -2.3863, -1.83553, -1.30092, -0.776683,
    -0.258268, 0.258268, 0.776683, 1.30092, 1.83553, 2.3863, 2.96138, 3.57377, 4.24812,
    5.04836}

For comparison purposes, let's compare that and a machine-precision evaluation of NSolve[] with an arbitrary-precision evaluation:

s2 = Sort[x /. NSolve[HermiteH[18, x], x]];

sb = Sort[x /. NSolve[HermiteH[18, x], x, WorkingPrecision -> 30]];

{Norm[sb - s1, ∞], Norm[sb - s2, ∞]}
   {4.44089*10^-15, 2.40696*10^-13}

and we see that the Jacobi-based method gives better results.

Zeros of high degree polynomials

3 Answers3