The interpolation problem - Fundamentals of Numerical Computation

In this chapter, we use $t_k$ for the nodes and $x$ to denote the continuous independent variable.

5.1.1Polynomials¶

Polynomials are the obvious first candidate to serve as interpolating functions. They are easy to work with, and in Polynomial interpolation we saw that a linear system of equations can be used to determine the coefficients of a polynomial that passes through every member of a set of given points in the plane. However, it’s not hard to find examples for which polynomial interpolation leads to unusable results.

Example 5.1.1 (Trouble in polynomial interpolation)

Julia

MATLAB

Python

Example 5.1.1

Here are some points that we could consider to be observations of an unknown function on $[-1,1]$ .

using Plots
n = 5
t = range(-1, 1, n+1)
y = @. t^2 + t + 0.05 * sin(20t)
scatter(t, y, label="data", legend=:top)

The polynomial interpolant, as computed using fit, looks very sensible. It’s the kind of function you’d take home to meet your parents.

using Polynomials
p = Polynomials.fit(t, y, n)     # interpolating polynomial
plot!(p, -1, 1, label="interpolant")

But now consider a different set of points generated in almost exactly the same way.

n = 18
t = range(-1, 1, n+1)
y = @. t^2 + t + 0.05 * sin(20t)
scatter(t, y, label="data", leg=:top)

The points themselves are unremarkable. But take a look at what happens to the polynomial interpolant.

p = Polynomials.fit(t, y, n)
x = range(-1, 1, 1000)    # use a lot of points
plot!(x, p.(x), label="interpolant")

Surely there must be functions that are more intuitively representative of those points!

Example 5.1.1

Here are some points that we could consider to be observations of an unknown function on $[-1,1]$ .

n = 5;
t = linspace(-1,1,n+1)';  
y = t.^2 + t + 0.05 * sin(20 * t);
clf, scatter(t,y)

The polynomial interpolant, as computed using polyfit, looks very sensible. It’s the kind of function you’d take home to meet your parents.

c = polyfit(t, y, n);     % polynomial coefficients
p = @(x) polyval(c, x);
hold on
fplot(p, [-1 1])
legend('data', 'interpolant', 'location', 'north');

But now consider a different set of points generated in almost exactly the same way.

n = 18;
t = linspace(-1, 1, n+1);
y = t.^2 + t + 0.05 * sin(20 * t);
clf, scatter(t, y)

The points themselves are unremarkable. But take a look at what happens to the polynomial interpolant.

c = polyfit(t, y, n);     % polynomial coefficients
p = @(x) polyval(c, x);
hold on, fplot(p, [-1 1])
legend('data', 'interpolant', 'location', 'north');

Surely there must be functions that are more intuitively representative of those points!

Example 5.1.1

Here are some points that we could consider to be observations of an unknown function on $[-1,1]$ .

n = 5
t = linspace(-1, 1, n + 1)
y = t**2 + t + 0.05 * sin(20 * t)
fig, ax = subplots()
plot(t, y, "o", label="data")
xlabel("$x$"),  ylabel("$y$");

The polynomial interpolant, as computed using fit, looks very sensible. It’s the kind of function you’d take home to meet your parents.

p = poly1d(polyfit(t, y, n))  # interpolating polynomial
tt = linspace(-1, 1, 400)
ax.plot(tt, p(tt), label="interpolant")
ax.legend()
fig

But now consider a different set of points generated in almost exactly the same way.

n = 18
t = linspace(-1, 1, n + 1)
y = t**2 + t + 0.05 * sin(20 * t)
fig, ax = subplots()
plot(t, y, "o", label="data")
xlabel("$x$"),  ylabel("$y$");

The points themselves are unremarkable. But take a look at what happens to the polynomial interpolant.

p = poly1d(polyfit(t, y, n))
ax.plot(tt, p(tt), label="interpolant")
ax.legend()
fig

Surely there must be functions that are more intuitively representative of those points!

In Chapter 9 we explore the large oscillations in the last figure of Example 5.1.1; it turns out that one must abandon either equally spaced nodes or $n\to\infty$ for polynomials. In the rest of this chapter we will keep $n$ fairly small and let the nodes be unrestricted.

5.1.2Piecewise polynomials¶

In order to keep polynomial degrees small while interpolating large data sets, we will choose interpolants from the piecewise polynomials. Specifically, the interpolant $p$ must be a polynomial on each subinterval $[t_{k-1},t_k]$ for $k=1,\ldots,n$ .

Usually we designate in advance a maximum degree $d$ for each polynomial piece of $p(x)$ . An important property of the piecewise polynomials of degree $d$ is that they form a vector space: that is, any linear combination of piecewise polynomials of degree $d$ is another piecewise polynomial of degree $d$ . If $p$ and $q$ share the same node set, then the combination is piecewise polynomial on that node set.

Example 5.1.3 (Piecewise polynomial interpolation)

Julia

MATLAB

Python

Example 5.1.3

Let us recall the data from Example 5.1.1.

n = 12
t = range(-1, 1, n+1)
y = @. t^2 + t + 0.5 * sin(20t)
scatter(t, y, label="data", leg=:top)

Here is an interpolant that is linear between each consecutive pair of nodes, using plinterp from Piecewise linear interpolation.

p = FNC.plinterp(t, y)
plot!(p, -1, 1, label="piecewise linear")

We may prefer a smoother interpolant that is piecewise cubic, generated using Spline1D from the Dierckx package.

using Dierckx
p = Spline1D(t, y)
plot!(x -> p(x), -1, 1, label="piecewise cubic")

Example 5.1.3

Let us recall the data from Example 5.1.1.

n = 18;
t = linspace(-1, 1, n+1);
y = t.^2 + t + 0.05 * sin(20 * t);
clf, scatter(t, y)

Here is an interpolant that is linear between each consecutive pair of nodes, using interp1 from MATLAB.

x = linspace(-1, 1, 400)';
hold on, plot(x, interp1(t, y, x))
title('Piecewise linear interpolant')

We may prefer a smoother interpolant that is piecewise cubic, generated using Spline1D from the Dierckx package.

cla
scatter(t, y)
plot(x, interp1(t, y, x, 'spline'))
title('Piecewise cubic interpolant')

Example 5.1.3

Let us recall the data from Example 5.1.1.

clf
n = 12
t = linspace(-1, 1, n + 1)
y = t**2 + t + 0.5 * sin(20 * t)
fig, ax = subplots()
scatter(t, y, label="data")
xlabel("$x$"),  ylabel("$y$");

Here is an interpolant that is linear between each consecutive pair of nodes, using plinterp from Piecewise linear interpolation.

from scipy.interpolate import interp1d
tt = linspace(-1, 1, 400)
p = interp1d(t, y, kind="linear")
ax.plot(tt, p(tt), label="piecewise linear")
ax.legend()
fig

We may prefer a smoother interpolant that is piecewise cubic:

scatter(t, y, label="data")
p = interp1d(t, y, kind="cubic")
tt = linspace(-1, 1, 400)
plot(tt, p(tt), label="cubic spline")
xlabel("$x$"),  ylabel("$y$")
legend();

We will consider piecewise linear interpolation in more detail in Piecewise linear interpolation, and we look at piecewise cubic interpolation in Cubic splines.

5.1.3Conditioning of interpolation¶

In the interpolation problem we are given the values $(t_k,y_k)$ for $k=0,\ldots,n$ . Let us consider the nodes $t_k$ of the problem to be fixed, and let $a=t_0$ , $b=t_n$ . Then the data for the interpolation problem consists of a vector $\mathbf{y}$ , and the result of the problem is a function on $[a,b]$ .

Let $\mathcal{I}$ be a prescription for producing the interpolant from a data vector. That is, $\mathcal{I}(\mathbf{y})=p$ , where $p(t_k)=y_k$ for all $k$ . The interpolation methods we will consider are all linear, in the sense that

\cI(\alpha\mathbf{y} + \beta\mathbf{z}) = \alpha \cI(\mathbf{y}) + \beta \cI(\mathbf{z})

(5.1.1)

for all vectors $\mathbf{y},\mathbf{z}$ and scalars $\alpha,\beta$ .

Linearity greatly simplifies the analysis of interpolation. To begin with, for any data vector $\mathbf{y}$ we have the standard expression $\mathbf{y}=\sum y_k \mathbf{e}_k$ , where, as always, $\mathbf{e}_k$ is a column of an identity matrix.^[1] Hence by linearity,

\cI( \mathbf{y} ) = \cI \left( \sum_{k=0}^n y_k \mathbf{e}_k \right) = \sum_{k=0}^n y_k \cI( \mathbf{e}_k ).

(5.1.2)

The functions appearing within the sum above have particular significance.

For any set of $n+1$ nodes, there are $n+1$ cardinal functions $\phi_0,\ldots,\phi_n$ , each singling out a different interpolation node in the set. We finish (5.1.2) by writing

\cI( \mathbf{y} ) = \sum_{k=0}^n y_k \phi_k.

(5.1.3)

In the following result we use the function infinity-norm or max-norm defined by

\| f\|_{\infty} = \max_{x \in [a,b]} |f(x)|.

(5.1.4)

Proof

Suppose the data vector is perturbed from $\mathbf{y}$ to $\mathbf{y}+ \mathbf{d}$ . Then

\cI(\mathbf{y} + \mathbf{d}) - \cI(\mathbf{y}) = \cI(\mathbf{d}) = \sum_{k=0}^n d_k \phi_k.

(5.1.6)

Hence

\frac{\bigl\|\cI(\mathbf{y} + \mathbf{d}) - \cI(\mathbf{y}) \bigr\|_{\infty}}{\| \mathbf{d} \|_{\infty}} = \left\|\, \sum_{k=0}^{n} \frac{d_k}{\|\mathbf{d} \|_{\infty}} \phi_k \, \right\|_{\infty}.

(5.1.7)

The absolute condition number maximizes this quantity over all $\mathbf{d}$ . Suppose $j$ is such that $\|\phi_j\|_\infty$ is maximal. Then let $\mathbf{d}=\mathbf{e}_j$ and the first inequality in (5.1.5) follows. The other inequality follows from the triangle inequality:

\left\| \, \sum_{k=0}^{n} \frac{d_k}{\|\mathbf{d} \|_{\infty}} \phi_k \, \right\|_{\infty} \le \sum_{k=0}^{n} \frac{|d_k|}{\|\mathbf{d} \|_{\infty}} \| \phi_k \|_\infty.

(5.1.8)

Since $|d_k|\le \|\mathbf{d}\|_\infty$ for all $k$ , this finishes (5.1.5).

Example 5.1.4 (Conditioning of interpolation)

Julia

MATLAB

Python

Example 5.1.4

In Example 5.1.3 we saw a big difference between polynomial interpolation and piecewise polynomial interpolation of some arbitrarily chosen data. The same effects can be seen clearly in the cardinal functions, which are closely tied to the condition numbers.

n = 18
t = range(-1, 1, n+1)
y = [zeros(9); 1; zeros(n - 9)];  # data for 10th cardinal function

scatter(t, y, label="data")

ϕ = Spline1D(t, y)
plot!(x -> ϕ(x), -1, 1;
    label="spline",
    xlabel=L"x",  ylabel=L"\phi(x)",
    title="Piecewise cubic cardinal function")

The piecewise cubic cardinal function is nowhere greater than one in absolute value. This happens to be true for all the cardinal functions, ensuring a good condition number for any interpolation with these functions. But the story for global polynomials is very different.

scatter(t, y, label="data")

ϕ = Polynomials.fit(t, y, n)
plot!(x -> ϕ(x), -1, 1;
    label="polynomial",  legend=:top,
    xlabel=L"x",  ylabel=L"\phi(x)", 
    title="Polynomial cardinal function")

From the figure we can see that the condition number for polynomial interpolation on these nodes is at least 500.

Example 5.1.4

n = 18;
t = linspace(-1, 1, n+1)';
y = [zeros(9, 1); 1; zeros(n - 9, 1)];    % 10th cardinal function
clf, scatter(t, y)
hold on
x = linspace(-1, 1, 400)';
plot(x, interp1(t, y, x, 'spline'))
title('Piecewise cubic cardinal function') 
xlabel('x'), ylabel('p(x)')

clf, scatter(t, y)
c = polyfit(t, y, n);
hold on, plot(x, polyval(c, x))
title('Polynomial cardinal function')
xlabel('x'), ylabel(('p(x)'));

From the figure we can see that the condition number for polynomial interpolation on these nodes is at least 500.

Example 5.1.4

clf
n = 18
t = linspace(-1, 1, n + 1)
y = zeros(n + 1)
y[9] = 1.0
p = interp1d(t, y, kind="cubic")

scatter(t, y, label="data")
tt = linspace(-1, 1, 400)
plot(tt, p(tt), label="cardinal function")
title("Cubic spline cardinal function")
legend();

p = poly1d(polyfit(t, y, n))
scatter(t, y, label="data")
plot(tt, p(tt), label="cardinal function")
xlabel("$x$")
ylabel("$y$")
title("Polynomial cardinal function")
legend();

From the figure we can see that the condition number for polynomial interpolation on these nodes is at least 500.

5.1.4Exercises¶

Exercise 5.1.1

⌨ Create data by entering

t = -2:4;  y = tanh.(t);

(a) Use fit to construct and plot the polynomial interpolant of the data, superimposed on a scatter plot of the data.

(b) Use Spline1D to construct and plot a piecewise cubic interpolant of the data, superimposed on a scatter plot of the data.

Exercise 5.1.2

⌨ The following table gives the life expectancy in the U.S. by year of birth.

1980	1985	1990	1995	2000	2005	2010
73.7	74.7	75.4	75.8	77.0	77.8	78.7

(a) Defining “year since 1980” as the independent variable, use fit to construct and plot the polynomial interpolant of the data.

(b) Use Spline1D to construct and plot a piecewise cubic interpolant of the data.

(c) Use both methods to estimate the life expectancy for a person born in 2007. Which value is more believable?

Exercise 5.1.3

⌨ The following two vectors define a flying saucer shape.

x = [ 0,0.51,0.96,1.06,1.29,1.55,1.73,2.13,2.61,
2.19,1.76,1.56,1.25,1.04,0.58,0 ]
y = [ 0,0.16,0.16,0.43,0.62,0.48,0.19,0.18,0,
-0.12,-0.12,-0.29,-0.30,-0.15,-0.16,0 ]

We can regard both $x$ and $y$ as functions of a parameter $s$ , with the points being values given at $s=0,1,\ldots,15$ .

(a) Use Spline1D once on each coordinate as functions of $s$ , and make a picture of the flying saucer.

(b) One drawback of the result in part (a) is the noticeable corner at the left side, which corresponds to $s=0$ from above and $s=15$ from below. There is a periodic variation on cubic spline interpolation that you can invoke by adding the keyword periodic=true to the Spline1D call. Use this to re-plot the flying saucer.

Footnotes¶

To be precise, we are using $\mathbf{e}_k$ to mean column number $k+1$ from an $(n+1)\times (n+1)$ identity matrix, since in linear algebra we start indexing at 1.
↩

Preface

5. Piecewise interpolation

Preface

Piecewise linear interpolation