Newton Polynomial Interpolation#

The method is similar to the Lagrange polynomial interpolation, just that the polynomial is computed differently. The end polynomial is similar, so it should give the same plot as the Lagrange polynomial.

\[f(x) = \sum _{i=0}^n a_i \, n_i (x)\]

where

\[n_i(x) = \prod _{j=0}^{i-1}(x - x_j)\]

The special feature of this polynomial is that the coefficents \(a_i\) can be easily determined.

\[f(x_0) = a_0 = y_0 f(x_1) = a_0 + a_1(x_1 - x_0) = y_1\]

which becomes

\[a_1 = \frac{y_1 - y_0}{x_1 - x_0}\]

now insert the point \((x_2,y_2)\) to find \(a_2\)

\[a_2 = \cfrac {\cfrac {y_2 - y_1}{x_2 - x_1} - \cfrac {y_1 - y_0}{x_1 - x_0}}{x_2 - x_0}\]

insert the point \((x_3,y_3)\) to find \(a_3\), produces another equation of divided differences.

\[\begin{split}f[x_1,x_0] &= \frac {y_1 - y_0}{x_1 - x_0} \\ f[x_2,x_1,x_0] &= \cfrac {\cfrac {y_2 - y_1}{x_2 - x_1} - \cfrac {y_1 - y_0}{x_1 - x_0}}{x_2 - x_0} \\ &= \frac {f[x_2,x_1] - f[x_1,x_0]}{x_2 - x_0}\end{split}\]

this produces the following iteration method

\[f[x_k,x_{k-1},...,x_1,x_0] = \frac {f[x_k,x_{k-1},...,x_2,x_1] - f[x_{k-1},x_{k-2},...,x_1,x_0]}{x_k - x_0}\]

produces the following divided difference coefficients in matrix form

\[\begin{split}\begin{array}{lcccc} &y_0 \; &f[x_1,x_0] \; &f[x_2,x_1,x_0] \; &f[x_3,x_2,x_1,x_0] \; &f[x_4,x_3,x_2,x_1,x_0] \\ &y_1 \; &f[x_2,x_1] \; &f[x_3,x_2,x_1] \; &f[x_4,x_3,x_2,x_1] \; &0\\ &y_2 \; &f[x_3,x_2] \; &f[x_4,x_3,x_2] \; &0 \; &0\\ &y_3 \; &f[x_4,x_3] \; &0 \; &0 \; &0\\ &y_4 \; &0 \; &0 \; &0 \; &0 \end{array}\end{split}\]

The first row gives all the polynomial coefficients \(a_0,a_1,a_2,a_3\)

../_images/newton_interp.png — Newton interpolation through 4 points#

The Newton and Lagrange interpolation are compared

# https://pythonnumericalmethods.berkeley.edu/notebooks/chapter17.05-Newtons-Polynomial-Interpolation.html
import numpy as np
import matplotlib.pyplot as plt
import seaborn as sns

sns.set_theme()

def divided_diff(x, y):
    '''
    function to calculate the divided
    differences table
    '''
    n = len(y)
    coef = np.zeros([n, n])
    # the first column is y
    coef[:,0] = y

    for j in range(1,n):
        for i in range(n-j):
            coef[i][j] = \
           (coef[i+1][j-1] - coef[i][j-1]) / (x[i+j]-x[i])

    return coef

def newton_poly(coef, x_data, x):
    '''
    evaluate the newton polynomial
    at x
    '''
    n = len(x_data) - 1
    p = coef[n]
    for k in range(1,n+1):
        p = coef[n-k] + (x -x_data[n-k])*p
    return p

class LagrangePoly:

    def __init__(self, X, Y):
        self.n = len(X)
        self.X = np.array(X)
        self.Y = np.array(Y)

    def basis(self, x, j):
        b = [(x - self.X[m]) / (self.X[j] - self.X[m])
             for m in range(self.n) if m != j]
        return np.prod(b, axis=0) * self.Y[j]

    def interpolate(self, x):
        b = [self.basis(x, j) for j in range(self.n)]
        return np.sum(b, axis=0)

x = np.array([-5, -1, 0, 2])
y = np.array([-2, 6, 1, 3])
# get the divided difference coef
a_s = divided_diff(x, y)[0, :]

# evaluate on new data points
x_new = np.linspace(max(x), min(x), int((max(x)-min(x))*10)) # -5, 2.1, .1
y_new = newton_poly(a_s, x, x_new)

lp = LagrangePoly(x, y)


plt.figure(figsize = (12, 8))
plt.plot(x, y, 'ro', label='data points')
plt.plot(x_new, y_new, 'y', label='Newton poly')
plt.plot(x_new, lp.interpolate(x_new), 'k', linestyle=':', label='Lagrange')
plt.legend()
plt.title('Newton Interpolation Polynomial\nCoincides with Lagrange Polynomial')
plt.xlabel('$x$')
plt.ylabel('$f(x)$')

plt.show()
#plt.savefig('../../figures/newton_interp.png')

As stated before both methods produce the same results. If we were to test it on more data points it could oscillate wildly at the ends (just like the lagrange polynomial), therefore the cubic spline method is shown next, which overcomes this problem.

Newton Polynomial Interpolation#

This Page