# Lagrange interpolation of 2000 equidistant points on $\frac{1}{1+x^2}$ crashes mathematica/does not evaluate

My Code:

Points := Table[-5 + 10 j/2048, {j, 1, 2047}]
Data := Table[1/(1 + x^2), {x, Points}]
Data2 := Transpose[{Points, Data}]
p[x_] = InterpolatingPolynomial[Data2, x]


I am not sure why this crashes Mathematica. It seems it cannot handle this. I have a HW assignment that requires for me to find the interpolating polynomial of $$\frac{1}{1+x^2}$$ at the points $$-5+10j/2048$$ for $$j$$ between 1 and 2047. The code above is my best attempt at doing so. I could not figure out how to make a table of the form $$\{\{x_1,y_1\}...\}$$ so I just combined the two sets by transposition. Mathematica handles The first three lines properly. However, when trying to evaluate

p[x_] = InterpolatingPolynomial[Data2, x]

after about 2-3 minutes crashes. I am not sure what the issue is. Is this impossible for mathematica? I tried the command parallelize but that di not help. The goal of the problem is to find the error of this interpolatoin, and compare it to the error of interpolation using chebyshev points (which I will try doing after I can do this)

• "... requires for me to find the interpolating polynomial..." -- I wonder what one would do with it after getting hold of it. In any case, try it with floating-point numbers, Points = Table[-5. + 10 j/2048, {j, 1, 2047}]. You'll run into a different problem, and then you can try it with arbitrary-precision floats, Points = Table[-5.16 + 10 j/2048, {j, 1, 2047}]. The numerators and denominators of the exact rational coefficients get very large. Oct 1 '20 at 20:29
• Oops, you'll need a lot more precision than 16 digits... Oct 1 '20 at 20:31
• @MichaelE2 I am to find its $\|\cdot\|_\infty$ error ands compare that error with interpolation but with chebyshev points. I tried the $-5.$ and it computed it! Could you explain what you did there? Why did it work now but not before? I tested a few values and it seems to work. It just tell me "General::stop: Further output of General::munfl will be suppressed during this calculation." as there are a lot of terms that don't really matter as they are so small. I think you solved mi issue though! Thank you Oct 1 '20 at 20:42
• @DanielHuber I think that's the point of the exercise. And to show that the use of the Chebyshev points + Chebyshev basis is stable to an extremely high degree and does not exhibit the bad behavior you describe. — I don't have time to answer, but someone should point out that InterpolatingPolynomial does Newton interpolation, not Lagrange. Oct 2 '20 at 16:38
• StatisticsLibraryBarycentricInterpolation gives a stable interpolation at a given set of nodes, esp. when the weights are known. Using it would show the issues with polynomial interpolation on a uniform grid (Runge phenomenon), without adding numerical error from evaluating the interpolating polynomial. You still need high precision. StatisticsLibraryBarycentricInterpolation[Points, Data, Weights -> N[(-1)^# Binomial[Length@Points - 1, #] &@ Range[0, Length@Points - 1], Max[Precision /@ Points]]][x]` -- See mathematica.stackexchange.com/a/191600/4999 for examples. Oct 2 '20 at 17:42