Polynomial interpolation - Fundamentals of Numerical Computation

The United States carries out a census of its population every 10 years. Suppose we want to know the population at times in between the census years, or to estimate future populations. One technique is to find a polynomial that passes through all the data points.^[1]

As posed in Definition 2.1.1, the polynomial interpolation problem has a unique solution. Once the interpolating polynomial is found, it can be evaluated anywhere to estimate or predict values.

2.1.1Interpolation as a linear system¶

Given data $(t_i,y_i)$ for $i=1,\ldots,n$ , we seek a polynomial

p(t) = c_1 + c_{2} t + c_3t^2 + \cdots + c_{n} t^{n-1},

(2.1.1)

such that $y_i=p(t_i)$ for all $i$ . These conditions are used to determine the coefficients $c_1\ldots,c_n$ :

\begin{split} c_1 + c_2 t_1 + \cdots + c_{n-1}t_1^{n-2} + c_nt_1^{n-1} &= y_1, \\ c_1 + c_2 t_2 + \cdots + c_{n-1}t_2^{n-2} + c_nt_2^{n-1} &= y_2, \\ c_1 + c_2 t_3 + \cdots + c_{n-1}t_3^{n-2} + c_nt_3^{n-1} &= y_3, \\ \vdots \qquad & \\ c_1 + c_2 t_n + \cdots + c_{n-1}t_n^{n-2} + c_nt_n^{n-1} &= y_n. \end{split}

(2.1.2)

These equations form a linear system for the coefficients $c_i$ :

\begin{bmatrix} 1 & t_1 & \cdots & t_1^{n-2} & t_1^{n-1} \\ 1 & t_2 & \cdots & t_2^{n-2} & t_2^{n-1} \\ 1 & t_3 & \cdots & t_3^{n-2} & t_3^{n-1} \\ \vdots & \vdots & & \vdots & \vdots \\ 1 & t_n & \cdots & t_n^{n-2} & t_n^{n-1} \\ \end{bmatrix} \begin{bmatrix} c_1 \\ c_2 \\ c_3 \\ \vdots \\ c_n \end{bmatrix} = \begin{bmatrix} y_1 \\ y_2 \\ y_3 \\ \vdots \\ y_n \end{bmatrix},

(2.1.3)

or more simply, $\mathbf{V} \mathbf{c} = \mathbf{y}$ . The matrix $\mathbf{V}$ is of a special type.

Polynomial interpolation can therefore be posed as a linear system of equations with a Vandermonde matrix.

Example 2.1.1 (Linear system for polynomial interpolation)

Julia

MATLAB

Python

Example 2.1.1

We create two vectors for data about the population of China. The first has the years of census data and the other has the population, in millions of people.

year = [1982, 2000, 2010, 2015]; 
pop = [1008.18, 1262.64, 1337.82, 1374.62];

It’s convenient to measure time in years since 1980. We use .- to subtract a scalar from every element of a vector. We will also use a floating-point value in the subtraction, so the result is also in double precision.

t = year .- 1980.0
y = pop;

Now we have four data points $(t_1,y_1),\dots,(t_4,y_4)$ , so $n=4$ and we seek an interpolating cubic polynomial. We construct the associated Vandermonde matrix:

V = [ t[i]^j for i in 1:4, j in 0:3 ]

4×4 Matrix{Float64}:
 1.0   2.0     4.0      8.0
 1.0  20.0   400.0   8000.0
 1.0  30.0   900.0  27000.0
 1.0  35.0  1225.0  42875.0

To solve for the vector of polynomial coefficients, we use a backslash to solve the linear system:

c = V \ y

4-element Vector{Float64}:
 962.2387878787875
  24.127754689754774
  -0.5922620490620537
   0.00684386724386731

The algorithms used by the backslash operator are the main topic of this chapter. As a check on the solution, we can compute the residual.

y - V * c

4-element Vector{Float64}:
 0.0
 0.0
 2.2737367544323206e-13
 0.0

Using floating-point arithmetic, it is not realistic to expect exact equality of quantities; a relative difference comparable to $\macheps$ is all we can look for.

By our definitions, the elements of c are coefficients in ascending-degree order for the interpolating polynomial. We can use the polynomial to estimate the population of China in 2005:

using Polynomials
p = Polynomial(c)    # construct a polynomial
p(2005-1980)         # include the 1980 time shift

1302.2043001443

The official population value for 2005 was 1303.72, so our result is rather good.

We can visualize the interpolation process. First, we plot the data as points.

using Plots
scatter(t, y;
    label="actual",  legend=:topleft,
    xlabel="years since 1980",  ylabel="population (millions)", 
    title="Population of China")

We want to superimpose a plot of the polynomial. We do that by evaluating it at a vector of points in the interval. The dot after the name of the polynomial is a universal way to apply a function to every element of an array, a technique known as broadcasting.

The range function constructs evenly spaced values given the endpoints and either the number of values, or the step size between them.

Adding a dot to the end of a function name causes it to be broadcast, i.e., applied to every element of a vector or matrix.

# Choose 500 times in the interval [0,35].
tt = range(0, 35, 500)

# Evaluate the polynomial at all the vector components.
yy = p.(tt)

foreach(println, yy[1:4])

962.2387878787875

963.9282039963299
965.6118068288089
967.2896105457506

Now we use plot! to add to the current plot, rather than replacing it.

plot!(tt, yy, label="interpolant")

Example 2.1.1

We create two column vectors for data about the population of China. The first has the years of census data and the other has the population, in millions of people.

year = [1982; 2000; 2010; 2015]; 
pop = [1008.18; 1262.64; 1337.82; 1374.62];

It’s convenient to measure time in years since 1980.

t = year - 1980;
y = pop;

Now we have four data points $(t_1,y_1),\dots,(t_4,y_4)$ , so $n=4$ and we seek an interpolating cubic polynomial. We construct the associated Vandermonde matrix:

V = vander(t)

To solve for the vector of polynomial coefficients, we use a backslash to solve the linear system:

c = V \ y

The algorithms used by the backslash operator are the main topic of this chapter. As a check on the solution, we can compute the residual.

y - V * c

Using floating-point arithmetic, it is not realistic to expect exact equality of quantities; a relative difference comparable to $\macheps$ is all we can look for.

By our definitions, the elements of c are coefficients in descending-degree order for the interpolating polynomial. We can use the polynomial to estimate the population of China in 2005:

p = @(t) polyval(c, t - 1980);  % include the 1980 time shift
p(2005)

The official population value for 2005 was 1303.72, so our result is rather good.

We can visualize the interpolation process. First, we plot the data as points.

scatter(year, y)
xlabel("years since 1980")
ylabel("population (millions)")
title("Population of China")

We want to superimpose a plot of the polynomial. We do that by evaluating it at a vector of points in the interval.

tt = linspace(1980, 2015, 500);    % 500 times in the interval [1980, 2015]
yy = p(tt);                        % evaluate p at all the vector elements
yy(1:4)

Now we use plot! to add to the current plot, rather than replacing it.

hold on 
plot(tt, yy)
legend("data", "interpolant", "location", "northwest");

Example 2.1.1

We create two vectors for data about the population of China. The first has the years of census data and the other has the population, in millions of people.

year = arange(1980, 2020, 10)   # from 1980 to 2020 by 10
pop = array([984.736, 1148.364, 1263.638, 1330.141])

It’s convenient to measure time in years since 1980.

t = year - 1980
y = pop

Now we have four data points $(t_1,y_1),\dots,(t_4,y_4)$ , so $n=4$ and we seek an interpolating cubic polynomial. We construct the associated Vandermonde matrix:

V = vander(t)
print(V)

[[    0     0     0     1]
 [ 1000   100    10     1]
 [ 8000   400    20     1]
 [27000   900    30     1]]

To solve a linear system $\mathbf{V} \mathbf{c} = \mathbf{y}$ for the vector of polynomial coefficients, we use solve (imported from numpy.linalg):

c = linalg.solve(V, y)
print(c)

[-6.95000e-05 -2.39685e-01  1.87666e+01  9.84736e+02]

The algorithms used by solve are the main topic of this chapter. As a check on the solution, we can compute the residual $\mathbf{y} - \mathbf{V} \mathbf{c}$ , which should be small (near machine precision).

print(y - V @ c)

[0. 0. 0. 0.]

By our definitions, the coefficients in c are given in descending order of power in $t$ . We can use the resulting polynomial to estimate the population of China in 2005:

p = poly1d(c)          # construct a polynomial
print(p(2005 - 1980))     # apply the 1980 time shift

1303.0119375

The official figure was 1303.72, so our result is rather good.

We can visualize the interpolation process. First, we plot the data as points. Then we add a plot of the interpolant, taking care to shift the $t$ variable back to actual years.

scatter(year, y, color="k", label="data");
tt = linspace(0, 30, 300)   # 300 times from 1980 to 2010
plot(1980 + tt, p(tt), label="interpolant");
xlabel("year");
ylabel("population (millions)");
title("Population of China");
legend();

2.1.2Exercises¶

Exercise 2.1.3

⌨ Here are population figures (in millions) for three countries over a 30-year period (from United Nations World Population Prospects, 2019).

	1990	2000	2010	2020
United States	252.120	281.711	309.011	331.003
India	873.278	1,056.576	1,234.281	1,380.004
Poland	37.960	38.557	38.330	37.847

(a) Use cubic polynomial interpolation to estimate the population of the USA in 2005.

(b) Use cubic polynomial interpolation to estimate when the population of Poland peaked during this time period.

(c) Use cubic polynomial interpolation to make a plot of the Indian population over this period. Your plot should be well labeled and show a smooth curve as well as the original data points.

Exercise 2.1.4

⌨ Here are the official population figures for the state of Delaware, USA, every ten years from 1790 to 1900: 59096, 64273, 72674, 72749, 76748, 78085, 91532, 112216, 125015, 146608, 168493, 184735. For this problem, use

t = \frac{\text{year} - 1860}{10}

(2.1.4)

as the independent (time) variable.

(a) Using only the data from years 1860 to 1900, plot the interpolating polynomial over the same range of years. Add the original data points to your plot as well.

(b) You might assume that adding more data will make the interpolation better. But this is not always the case. Use all the data above to create an interpolating polynomial of degree 11, and then plot that polynomial over the range 1860 to 1900. In what way is this fit clearly inferior to the previous one? (This phenomenon is studied in Chapter 9.)

Footnotes¶

We’re quite certain that the U.S. Census Bureau uses more sophisticated modeling techniques than the one we present here!
↩

Preface

2. Linear systems of equations

Preface

Computing with matrices