interpolation_details.ipynb

Kernel: Python 3 (ipykernel)

Appendix

Thecnical details of interpolation functions: interpolation_details.ipynb

Bibligraphy

• https://www.math.ust.hk/~mamu/courses/231/Slides/CH03_3A.pdf

• Narváez, Diana Marcela Devia. "Comparison between some techniques of interpolators: An application in engineering." Scientia et technica 24.1 (2019): 173-178. DOI: https://doi.org/10.22517/23447214.21341 (Open Access)

• Atkinson, Kendall E. An introduction to numerical analysis. John wiley & sons, 2008. PDF

Divided Differences

In spite of the good precision achieved by the Lagrange interpolating polynomials, analytical manipulation of such an expressions is rather complicated. Furthermore, when applying other polynomials-based techniques like Hermite polynomials, the algorithms present very different ways to achieve the final interpolation, making a comparison unclear.

Divided differences is a way to standardize the notation for interpolating polynomials. Suppose a polynomial $P_n(x)$ and write it in the next form:

$$P_n(x) = a_0 + a_1(x-x_0)+ a_2 (x-x_0)(x-x_1)+\cdots + a_n(x-x_0)\cdots (x-x_{n-1})$$

where $a_i$ are a set of constants to be determined from the given data $(x_i, y_i)$.

• Obtain $a_0$
  Note that due to the definition of an interpolant function, previous expression should satisfy:

$$P_{n}\left(x_{0}\right)=a_{0}=f\left(x_{0}\right)=y_{0}$$

• Obtain $a_1$ $$\frac{P_{n}(x)-a_{0}}{x-x_{0}}=\frac{P_{n}(x)-y_{0}}{x-x_{0}}=a_{1}+a_{2}\left(x-x_{1}\right)+\cdots+a_{n}\left(x-x_{1}\right) \cdots\left(x-x_{n}\right)$$

$$\frac{P_{n}\left(x_{1}\right)-y_{0}}{x_1-x_{0}}=a_{1}=\frac{y_{1}-y_{0}}{x_1-x_{0}}$$

• Obtain $a_2$ $$\frac{1}{x-x_{1}}\left(\frac{P_{n}(x)-y_{0}}{x-x_{0}}-\frac{y_{1}-y_{0}}{x_{1}-x_{0}}\right)=a_{2}+a_{3}\left(x-x_{2}\right)+\cdots+a_{n}\left(x-x_{2}\right) \cdots\left(x-x_{n}\right)$$ evaluating in $x=x_2$ $$a_{2}=\frac{1}{x_{2}-x_{1}}\left(\frac{y_{2}-y_{0}}{x_{2}-x_{0}}-\frac{y_{1}-y_{0}}{x_{1}-x_{0}}\right),$$ that we can rewrite like $$a_{2}=\frac{\frac{y_{2}-y_{1}}{x_{2}-x_{1}}-\frac{y_{1}-y_{0}}{x_{1}-x_{0}}}{x_{2}-x_{0}}.$$

In summary $$\begin{align} a_{0}=& f\left[x_{0}\right] \equiv y_{0} \\ a_{1}=& f\left[x_{0}, x_{1}\right] \equiv \frac{y_{1}-y_{0}}{x_{1}-x_{0}}\\ a_{2}=&f\left[x_{0}, x_{1}, x_{2}\right] \equiv \frac{\frac{y_{2}-y_{1}}{x_{2}-x_{1}}-\frac{y_{1}-y_{0}}{x_{1}-x_{0}}}{x_{2}-x_{0}} \end{align}$$

The key observation, by using the symmetry $f\left[x_{0}, x_{1}\right]=f\left[x_{1}, x_{0}\right]$, is $$f\left[x_{0}, x_{1}, x_{2}\right]=\frac{f\left[x_{1}, x_{2}\right]-f\left[x_{0}, x_{1}\right]}{x_{2}-x_{0}}$$

Allow us to definine the zeroth divided difference, $k=0$, of $x_i$ like

$$D_0[x_i] \equiv f[x_i] \equiv f(x_i) = y_i$$

the first divided difference, $k=1$ of $x_i$ like

$$D_1[x_i] \equiv f[x_i, x_{i+1}] \equiv \frac{f[x_{i+1}]-f[x_i]}{x_{i+1}-x_i}$$$$D_1[x_i] = \frac{D_{0}[x_{i+1}]-D_{0}[x_{i}]}{x_{i+1}-x_i}$$

the second divided difference, $k=2$ of $x_i$ like

$$D_2[x_i]\equiv f[x_i,x_{i+1},x_{i+2}] \equiv \frac{f[x_{i+1},x_{i+2}]-f[x_i,x_{i+1}]}{x_{i+2}-x_{i}}$$$$D_2[x_i]=\frac{D_1[x_{i+1}]-D_1[x_i]}{x_{i+2}-x_{i}}$$

successively until the kth divided difference

$$D_k[x_i] \equiv f[x_i, x_{i+1},\cdots, x_{i+k-1},x_{i+k}] \equiv \frac{f[x_{i+1},x_{i+2}\cdots, x_{i+k}]-f[x_i, x_{i+1},\cdots, x_{i+k-1}]}{x_{i+k}-x_i}$$$$D_k[x_i] = \frac{D_{k-1}[x_{i+1}]-D_{k-1}[x_{i}]}{x_{i+k}-x_i}$$

These expressions are the fundamental bricks for any interpolating method.

#Construction of a kth divided difference (recursive code)
def D( i, k, Xn, Yn ):
    #If k+i>N
    if i+k>=len(Xn):
        return 0
    #Zeroth divided difference
    elif k == 0:
        return Yn[i]
    #If higher divided difference
    else:
        return (D(i+1, k-1, Xn, Yn)-D(i, k-1, Xn, Yn))/(Xn[i+k]-Xn[i])

Example 2

As an example, Lagrange interpolation can be also derived by using divided differences, which is reached through the next equation:

$$P_n(x) = D_0[x_0] + \sum_{k=1}^n D_k[x_0] (x-x_0) \cdots (x-x_{k-1})$$

Note this expression is by far easier to be manipulated analytically as we can know the coefficients of each order.

Activity

Using the previous expression and the defined function for divided differences, show both methods to calculate Lagrange interpolators are equivalents.

---

Hermite Interpolation

From calculus we know that Taylor polynomials expand a function at a specific point $x_i$, being both functions (the original one and the Taylor function) exactly equal at any derivative-order at that point. Also, as mentioned before, a Lagrange polynomial, given a set of data points, passes through all those points at the same time. However if those points come from an unknown underlying function $f(x)$, the interpolant polynomial might (surely) differ from the real function at any superior derivative-order. So we have:

• Taylor polynomials are exact at any order, but that only remains true at a specific point.

• Lagrange polynomials pass through all points of a give dataset, but only at zeroth-order. Derivatives are not longer equal.

Once established these differences, we can introduce Hermite polynomials just as a generalization of both, Taylor and Lagrange polynomials.

At first, Hermite polynomials can be approximated at any desired order at all the points, as long as one has all these information. However, for the sake of simplicity and without loss of generality, we shall assume Hermite polynomials equal to the real function at zeroth and first-derivative order. For that, we also need to know the points associated to the first-derivate.

Let's suppose a dataset $\{x_i\}_i$ for $i = 0,1,\cdots,n$ with the respective values $\{f(x_i)\}_i$ and $\{f'(x_i)\}_i$. If we assume two different polynomials to fit each set of data, i.e. a polynomial for $\{f(x_i)\}_i$ and another for $\{f'(x_i)\}_i$, we obtain $2n+2$ coefficients, however zeroth-order coefficients can be put together so finally there are $2n+1$ independet coefficients to be determined. In this case, we assign the respective Hermite polynimial as $H_{2n+1}(x)$.

Derivation in terms of divided differences

Remembering the divided differences expression for a Lagrange polynomial

$$P_n(x) = D_0[x_0] + \sum_{k=1}^n D_k[x_0] (x-x_0) \cdots (x-x_{k-1})$$

and by defining a new sequence $\{z_0, z_1, \cdots, z_{2n+1}\}$ such that

For even $i$ $$z_{2i} = x_i \mbox{ for } i = 0,2,\cdots$$ while for odd $i$ $$z_{2i+1} = x_i \mbox{ for } i = 1,3,\cdots$$

However, divided differences has to be modified in order to include first-order derivatives:

Note that $f[z_0,z_1]$ sould be originally

$$f[z_0,z_1] = \frac{f[z_1]-f[z_0]}{z_1-z_0}$$

but replacing $z_0 = z_1 = x_0$ this would lead an indetermination. In order to solve this issue, this indertemination can be readily approximated to the derivative at $x_0=z_0$, so

$$f[z_0,z_1] = f'(x_0)$$

or using the previously defined notation

$$D_1[z_0] = f'(x_0)$$

Generally, for first-order divided differences we will have

even $i$ $$D_1[z_{2i}] = f'(x_i)$$ odd $i$ $$D_1[z_{2i+1}] = f'(x_i)$$

Higher order divided differences are calculated as usual: $$D_k[z_i] = \frac{D_{k-1}[z_{i+1}]-D_{k-1}[z_{i}]}{z_{i+k}-z_i}$$

Finally, the Hermite polynomial is built using the next expression

• Order zero coefficient: $$D_0[z_0]$$

• Hermite polynomials $$H_{2n+1}(x) = D_0[z_0] + \sum_{k=1}^{2n+1} D_k[z_0] (x-z_0) \cdots (x-z_{k-1})$$

Example

%pylab inline

Populating the interactive namespace from numpy and matplotlib

/usr/local/lib/python3.9/dist-packages/IPython/core/magics/pylab.py:159: UserWarning: pylab import has clobbered these variables: ['f']
`%matplotlib` prevents importing * from pylab and numpy
  warn("pylab import has clobbered these variables: %s"  % clobbered +

def f(x):
    '''
    WARNING: all the parts of the function must be 
    numpy arrays in order to get the element by element sum.
    '''
    x=np.asarray(x) # Force x → array
    return np.exp(2*x)-x

def fp(x):
    '''
    Derivative of f(x)
    '''
    x=np.asarray(x) # Force x → array
    return 2*np.exp(2*x)-1

from scipy import interpolate

x=[1,1.25,1.6]

H=interpolate.CubicHermiteSpline(x,f(x),fp(x))

X=np.linspace( x[0],x[-1]  )
plt.plot(X,f(X),'k-',lw=1,label="$f(x)$")
plt.plot(x,f(x),'r*')
plt.plot(X,H(X),'c:',lw=3,label='$H(x)$')

plt.plot(X,fp(X),'r-',lw=1,label="$f'(x)$")
plt.plot(x,fp(x),'r*')
plt.plot(X,H.derivative()(X),'g:',lw=3,label="$H'(x)$")


plt.grid()
plt.xlabel('$x$',size=12)
plt.legend()

<matplotlib.legend.Legend at 0x7f4af01582b0>

Example 3

Define a routine to calculate divided differences for Hermite polynomials.

#Construction of a kth divided difference for Hermite polynomials (recursive code)
def Dh( j, k, Zn, Yn, Ypn ):
    #If k+j>N
    if j+k>=len(Zn):
        return 0
    #Zeroth divided difference
    elif k == 0:
        return Yn[j/2]
    #First order divided difference (even indexes)
    elif k == 1 and j%2 == 0:
        return Ypn[j/2]
    #If higher divided difference
    else:
        return (Dh(j+1, k-1, Zn, Yn, Ypn)-Dh(j, k-1, Zn, Yn, Ypn))/(Zn[j+k]-Zn[j])

Activity

Calculate a routine, using the previous program for divided differences, that computes the Hermite polynomial given a dataset.

Generate a set of $N$ points of the function $\sin^2(x)$ between $0$ and $2\pi$, including an array of $x$ positions, $y = f(x)$ and first derivative $y' = f'(x)$.

Show which polynomial gives the best approximation to the real function, Hermite or Lagrange polynomial.

Solution:

nbviewer.ipython.org/github/sbustamante/ComputationalMethods/blob/master/activities/hermite-and-lagrange.ipynb

---