# Get the roots of a Hermite interpolating polynomial

I am using Python 3.7 to write a program that requires me to calculate the root of the Hermite interpolating polynomial, given two points $$\epsilon_0$$, $$\epsilon_1$$, the function ($$d(\epsilon_0)$$,$$d(\epsilon_1)$$) and the derivative values ($$d'(\epsilon_1)$$, $$d'(\epsilon_1)$$) at those points. I am using Scipy v1.3.0 and using the CubicHermiteSpline function from the scipy.interpolate library. The relevant extracts from the code are:

import numpy as np
from scipy.interpolate import BPoly,CubicHermiteSpline

#somewhere below inside a while loop with a counter variable k is this part

r=CubicHermiteSpline(eps[k-1:k+1],abs(l[k-1:k+1]), d1[k-1:k+1]).roots()
epsk=(np.abs(r - eps[k])).argmin()


Whereabs(l) contains the values for the polynomial and d1 contains the derivative values. The problem is the .roots() returns an empty array for the interval ($\epsilon_0$,$\epsilon_1$).

ValueError: attempt to get argmin of an empty sequence


This is because the interpolated polynomial from this interval looks like this:

How do I get all the three roots of the interpolated polynomial, which may not necessarily be inside the interval?

Edit: The numerical values: $$d(\epsilon_0)=1.00000188\\ d(\epsilon_1)=1.09393556\\ d'(\epsilon_0)=-4.30116854\\ d'(\epsilon_1)=-4.30428889$$ Find the roots of the hermite intpolation polynomial. Interpolated polynomial graph:

• Did you try creating an CubicSpline interpolator first and then using the roots method from there? – nicoguaro Jul 12 '19 at 17:02
• Should I? I can't find a way to introduce the derivatives. BPoly.from_derivatives has though, but a different object with no root-finding methods attached. – mm-crj Jul 12 '19 at 17:09
• It's a workaround. You first is the derivatives and get a CubicHermiteSpline and then you evaluate that in four points inside the domain and get a CubicSpline with that. – nicoguaro Jul 12 '19 at 17:13
• And the cubicspline will give the desired roots? Its the same Ppoly object though. Let me try.. :) – mm-crj Jul 12 '19 at 17:16

The interpolated polynomial does not have roots. Considering that the behavior outside the interpolation region holds is termed extrapolation.

You can explicitly use the polynomial, given by (as I explained in this post)

$$f(x) \approx N_1(x) u_1 + N_2(x) u_2 + |J|(N_3(x) u'_1 + N_4(x) u'_2)\quad \forall x\in [a, b]$$

with $$|J| = (b - a)/2$$ the Jacobian determinant of the transformation, and, for $$a=-1, b=1$$,

\begin{align} N_1 (x) &= \frac{1}{4} (x - 1)^2 (2 + x)\\ N_2 (x) &= \frac{1}{4} (x + 1)^2 (2 - x)\\ N_3 (x) &= \frac{1}{4} (x - 1)^2 (x + 1)\\ N_4 (x) &= \frac{1}{4} (x + 1)^2 (x - 1)\, . \end{align}

Then, you can use the general formula for cubic equations or use a method like Newton-Raphson.

• I saw it before, but it seemed like too much for such a simple task. Because in theory I already have the polynomial, just need to remove the endpoints. Thanks though, your effort is appreciated. – mm-crj Jul 12 '19 at 16:06
• @mm-crj, I don't have version 1.3.1 of SciPy so I don't have CubicHermiteSpline. Nevertheless, I tried with CubicSpline and the extrapolation is set to True by default. – nicoguaro Jul 12 '19 at 16:32
• Yeah, in the documentation for CubicHermiteSpline extrapolation is set to True by default and its giving values outside the range too. And I have found the roots by plotting it. But .roots or . solve returns an empty array. – mm-crj Jul 12 '19 at 16:48
• I also didn't want to use CubicHermiteSpline because of this version problem. but I found no other methods taking the derivatives into account(Bpoly.from derivatives) and also has a roots or solution method with it. – mm-crj Jul 12 '19 at 16:51