Polynomial Interpolation

5. Polynomial Interpolation#

In this chapter we provide basic implementation of polynomial interpolation methods.

5.1. Polynomials in Newton’s form#

We begin by providing the implementation of the nested-multiplication method for the evaluation of a polynomial at a given input point.

def evaluate_by_nested_multiplication(t:float,x:list,c:list):
    '''
    Evaluates the polynomial at a given point t.
    :param t: the point at which the polynomial must be evaluated.
    :param x: the n+1-dimensional vector of x, [x_0,\ldots,x_n]
    :param c: the n+1-dimensional vector of coefficients, [c_0,\ldots,c_n]
    :return: the value of the polynomial at t.
    '''
    # Identifies the lenght of x
    n = len(x)  # n+1 dimensional

    # Initializes u=c_n
    u = c[n - 1]
    # For i = n-1,...,0 (with decrements of -1)
    for i in range((n - 1) - 1, 0 - 1, -1):
        u = c[i] + (t - x[i]) * u
    return u

5.1.1. Example 1#

Assume we have a polynomial of degree \(n=2\) in Newton’s form, defined by the coefficients \(c=[3.75,-1.5,1]\) and by the values \(x=[-1/2,0,1/2,2]\). The polynomial is

\[p_2(x)=3.75-1.5(x+0.5)+(x-0)(x+0.5)\]

Let us evaluate the polynomial at the \(x\) points.

X = [-1/2,0,1/2]
c = [3.75,-1.5,1]
for t in X:
    print(f'p_2({t})={evaluate_by_nested_multiplication(t,X,c)}')

p_2(-0.5)=3.75
p_2(0)=3.0
p_2(0.5)=2.75

5.2. Computation of the coefficients by nested multiplication#

We now turn our attention to the computation of the coefficients. We extend the nested multiplication algorithm to compute the coefficients.

def coefficients_by_nested_multiplication(x:list, y:list):
    '''
    Computes the coefficients of the Newton interpolating polynomial.

    :param x: List of x-coordinates of the data points.
    :param y: List of y-coordinates of the data points.
    :return: List of coefficients of the Newton interpolating polynomial.
    '''

    # Identifies the number of points
    n = len(x)
    # Initializes the vector of coefficients
    c = [0] * n
    # Inizializes c_0
    c[0] = y[0]
    # For k = 1,...,n
    for k in range(1, n):
        # Populates d
        d = x[k] - x[k - 1]
        # Populates u
        u = c[k - 1]
        # For i = k-2,...,0
        for i in range(k - 2, -1, -1):
            a = x[k] - x[i]
            u = c[i] + u * a
            d = a * d
        c[k] = (y[k] - u) / d

    return c

5.2.1. Example 1#

Let us use the same example we have seen for the evaluation of a given polynomial. This time, assume we are give the points \(x=[-0.5, 0, 0.5, 2]\) and \(y=[3.75, 3, 2.75, 5]\) and we wish to find the interpolating polynomial in Newton’s form.

x = [-0.5, 0, 0.5, 2] 
y = [3.75, 3, 2.75, 5]

# Let us find the coefficients of the Newton's polynomial which interpolates the points.
c = coefficients_by_nested_multiplication(x,y)
print(c)

[3.75, -1.5, 1.0, 0.0]

Let us now verify that the Newton’s polynomial defined by these coefficients indeed interpolates the three points.

for t in x:
  print(f'p({t})={evaluate_by_nested_multiplication(t,x,c)}')

p(-0.5)=3.75
p(0)=3.0
p(0.5)=2.75
p(2)=5.0

5.3. Computation of the coefficients by the divided differences method#

We proceed by implementing the divided differences method.

import numpy as np
def coefficients_by_divided_differences(x:list,y:list):
    '''
    Computes the coefficients of the Newton interpolating polynomial using the divided differences table.

    :param x: List of x-coordinates of the data points.
    :param y: List of y-coordinates of the data points.
    :return: List of coefficients of the Newton interpolating polynomial.
    '''

    # Finds the number of interpolation points
    n = len(x)

    # Initializes the divided differences table as an np.array of zeros.
    f = np.zeros((n, n))

    # Sets the first column of the table
    for row in range(0,n):
        f[row,0] = y[row]

    # Populates the remaining columns
    for column in range(1,n):
        for row in range(0,n-column):
            f[row,column] = (f[row+1,column-1] - f[row,column-1])/(x[column+row]-x[row])
    # Returns the divided differences table.
    return f

Observe that the method returns the entire divided differences table. This can be useful in case we later want to extend the polynomial to interpolate an additional point.

5.3.1. Example 1#

Let us consider the same example we used with the nested multiplication method. This time we apply the divided differences method.

x = [-0.5, 0, 0.5, 2] 
y = [3.75, 3, 2.75, 5]

# Let us find the coefficients of the Newton's polynomial which interpolates the points.
c = coefficients_by_divided_differences(x,y)
print(c)

[[ 3.75 -1.5   1.    0.  ]
 [ 3.   -0.5   1.    0.  ]
 [ 2.75  1.5   0.    0.  ]
 [ 5.    0.    0.    0.  ]]

The coefficients of the polynomial are given in the first row of the divided differences table.

c[0]

array([ 3.75, -1.5 ,  1.  ,  0.  ])

5.4. Exercises#

Exercise 5.1

Assume you have found the divided differences table for the polynomial which interpolates the points with coordinates \(x=[x_0,\ldots,x_n]\), \(Y=[y_0,\ldots,y_n]\). Implement a function that takes the following input

  • The divided differences

  • The vector \(x\)

  • An additional point \((x_{n+1},y_{n+1})\). The function must the divided differences table for the polynomial which interpolates the points with coordinates \([x_0,\ldots,x_n,x_{n+1}]\), \(Y=[y_0,\ldots,y_n,y_{n+1}]\).

Contents