# Zeros of high degree polynomials

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 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 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 at 4:48 • Possible duplicate: mathematica.stackexchange.com/q/51098/4999 – Michael E2 Jun 29 at 4:50 ## 2 Answers 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 at 16:08

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