3

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

Mainak
  • 181
  • 5
  • 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. – Mainak 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.. :) – Mainak Jul 12 '19 at 17:16

1 Answers1

6

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.

nicoguaro
  • 8,500
  • 6
  • 23
  • 49
  • 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. – Mainak Jul 12 '19 at 16:06
  • 2
    @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. – Mainak 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. – Mainak Jul 12 '19 at 16:51