Lecture notes for TMA4125/4130/4135 Mathematics 4N/D

Polynomial interpolation: Newton interpolation

Anne Kværnø (modified by André Massing)

Jan 18, 2021

The Python codes for this note are given in polynomialinterpolation.py.

1 Introduction
We continue with an alternative approach to find the interpolation polynomial.

1.1 Newton interpolation

This is an alternative approach to find the interpolation polynomial. Let x0 , x1 , . . . , xn be n + 1 distinct
real numbers. Instead of using the Lagrange polynomials to write the interpolation polynomial in Lagrange
form, we will now employ the Newton polynomials ωi , i = 0, . . . , n. The Newton polynomials are defined
by as follows:

ω0 (x) = 1,
ω1 (x) = (x − x0 ),
ω2 (x) = (x − x0 )(x − x1 ),
...
ωn (x) = (x − x0 )(x − x1 ) · · · (x − xn−1 ),

or in more compact notation

i−1
Y
ωi (x) = (x − xk ). (1)
k=0

The so-called Newton form of a polynomial of degree n is an expansion of the form

n
X
p(x) = ci ωi (x)
i=0

or more explicitly

p(x) = cn (x − x0 )(x − x1 ) · · · (x − xn−1 ) + cn−1 (x − x0 )(x − x1 ) · · · (x − xn−2 ) + · · · + c1 (x − x0 ) + c0 .

In the light of this form of writing a polynomial, the polynomial interpolation problem leads to the
following observations. Let us start with a single node x0 , then f (x0 ) = p(x0 ) = c0 . Going one step further
and consider two nodes x0 , x1 . Then we see that f (x0 ) = p(x0 ) = c0 and f (x1 ) = p(x1 ) = c0 + c1 (x1 − x0 ).
The latter implies that the coefficient
f (x1 ) − f (x0 )
c1 = .
x1 − x0
Given three nodes x0 , x1 , x2 yields the coefficients c0 , c1 as defined above, and from

f (x2 ) = p(x2 ) = c0 + c1 (x2 − x0 ) + c2 (x2 − x0 )(x2 − x1 )

we deduce the coefficient

f (x1 )−f (x0 )

f (x2 ) − f (x0 ) − x1 −x0 (x2 − x0 )
c2 = .
(x2 − x0 )(x2 − x1 )

Playing with this quotient gives the much more structured expression

f (x2 )−f (x1 ) f (x1 )−f (x0 )

x2 −x1 − x1 −x0
c2 = .
(x2 − x0 )

This procedure can be continued and yields a so-called triangular systems that permits to define the
remaining coefficients c3 , . . . , cn . One sees quickly that the coefficient ck only depends on the interpolation
points (x0 , y0 ), . . . , (xk , yk ), where yi := f (xi ), i = 0, . . . , n.
We introduce the folllwing so-called finite difference notation for a function f . The 0th order finite
difference is defined to be f [x0 ] := f (x0 ). The 1st order finite difference is

f (x1 ) − f (x0 )
f [x0 , x1 ] := .
x1 − x0

The second order finite difference is defined by

f [x1 , x2 ] − f [x0 , x1 ]
f [x0 , x1 , x2 ] := .
x2 − x0

In general, the nth order finite difference of the function f , also called the nth Newton divided
difference, is defined recursively by
f [x1 , . . . , xn ] − f [x0 , . . . , xn−1 ]
f [x0 , . . . , xn ] := .
xn − x0

Newton’s method to solve the polynomial interpolation problem can be summarized as follows. Given
n + 1 interpolation points (x0 , y0 ), . . . , (xn , yn ), yi := f (xi ). If the order n interpolation polynomial is
expressed in Newton’s form

pn (x) = cn (x − x0 )(x − x1 ) · · · (x − xn−1 ) + cn−1 (x − x0 )(x − x1 ) · · · (x − xn−2 ) + · · · + c1 (x − x0 ) + c0 ,

then the coefficients

ck = f [x0 , . . . , xk ]

for k = 0, . . . , n. In fact, a recursion is in place

pn (x) = pn−1 (x) + f [x0 , . . . , xn ](x − x0 )(x − x1 ) · · · (x − xn−1 )

It is common to write the finite differences in a table, which for n = 3 will look like:

x0 f [x0 ]
f [x0 , x1 ]
x1 f [x1 ] f [x0 , x1 , x2 ]
f [x1 , x2 ] f [x0 , x1 , x2 , x3 ]
x2 f [x2 ] f [x1 , x2 , x3 ]
f [x2 , x3 ]
x3 f [x3 ]

2
Example 1 again: Given the points in Example 1. The corresponding table of divided differences
becomes:
0 1
−3/4
2/3 1/2 −3/4
−3/2
1 0

and the interpolation polynomial becomes

3 3 2 1 3
p2 (x) = 1 − (x − 0) − (x − 0)(x − ) = 1 − x − x2 .
4 4 3 4 4

1.2 Implementation
The method above is implemented as two functions:
• divdiff(xdata, ydata): Create the table of divided differences
• newtonInterpolation(F, xdata, x): Evaluate the interpolation polynomial.
Here, xdata and ydata are arrays with the interpolation points, and x is an array of values in which the
polynomial is evaluated.
def divdiff(xdata,ydata):
# Create the table of divided differences based
# on the data in the arrays x_data and y_data.
n = len(xdata)
F = np.zeros((n,n))
F[:,0] = ydata # Array for the divided differences
for j in range(n):
for i in range(n-j-1):
F[i,j+1] = (F[i+1,j]-F[i,j])/(xdata[i+j+1]-xdata[i])
return F # Return all of F for inspection.
# Only the first row is necessary for the
# polynomial.

def newton_interpolation(F, xdata, x):

# The Newton interpolation polynomial evaluated in x.
n, m = np.shape(F)
xpoly = np.ones(len(x)) # (x-x[0])(x-x[1])...
newton_poly = F[0,0]*np.ones(len(x)) # The Newton polynomial
for j in range(n-1):
xpoly = xpoly*(x-xdata[j])
newton_poly = newton_poly + F[0,j+1]*xpoly
return newton_poly

Run the code on the example above:

# Example: Use of divided differences and the Newton interpolation
# formula.
xdata = [0, 2/3, 1]
ydata = [1, 1/2, 0]
F = divdiff(xdata, ydata) # The table of divided differences
print(’The table of divided differences:\n’,F)

x = np.linspace(0, 1, 101) # The x-values in which the polynomial is evaluated

p = newton_interpolation(F, xdata, x)
plt.plot(x, p) # Plot the polynomial
plt.plot(xdata, ydata, ’o’) # Plot the interpolation points
plt.title(’The interpolation polynomial p(x)’)
plt.grid(True)
plt.xlabel(’x’);