3
$\begingroup$

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: enter image description here

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: enter image description here

$\endgroup$
  • $\begingroup$ Did you try creating an CubicSpline interpolator first and then using the roots method from there? $\endgroup$ – nicoguaro Jul 12 at 17:02
  • $\begingroup$ 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. $\endgroup$ – mm-crj Jul 12 at 17:09
  • $\begingroup$ 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. $\endgroup$ – nicoguaro Jul 12 at 17:13
  • $\begingroup$ And the cubicspline will give the desired roots? Its the same Ppoly object though. Let me try.. :) $\endgroup$ – mm-crj Jul 12 at 17:16
6
$\begingroup$

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.

$\endgroup$
  • $\begingroup$ 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. $\endgroup$ – mm-crj Jul 12 at 16:06
  • 2
    $\begingroup$ @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. $\endgroup$ – nicoguaro Jul 12 at 16:32
  • $\begingroup$ 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. $\endgroup$ – mm-crj Jul 12 at 16:48
  • $\begingroup$ 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. $\endgroup$ – mm-crj Jul 12 at 16:51

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.