Piecewise linear interpolation - Fundamentals of Numerical Computation

Piecewise linear interpolation is simply a game of connect-the-dots. That is, the data points are joined pairwise by line segments.

It should be clear from (5.2.1) that on each interval $[t_k,t_{k+1}]$ , $p(x)$ is a linear function passing through both $(t_k,y_k)$ and $(t_{k+1},y_{k+1})$ .

5.2.1Hat functions¶

Rather than basing an implementation on (5.2.1), we return to the idea used in Example 2.1.1 of choosing the interpolant from among the linear combinations of a preselected finite set of functions. In the present context we use, for $k=0,\ldots,n$ ,

H_k(x) = \begin{cases} \dfrac{x-t_{k-1}}{t_k-t_{k-1}} & \text{if $x\in[t_{k-1},t_k]$},\\[2.5ex] \dfrac{t_{k+1}-x}{t_{k+1}-t_{k}} & \text{if $x\in[t_{k},t_{k+1}]$},\\[2.5ex] 0 & \text{otherwise}. \end{cases} \qquad

(5.2.2)

The functions $H_0,\ldots,H_n$ are called hat functions. They depend on the node vector $\mathbf{t}$ , but this dependence is not usually indicated explicitly.

Each hat function is globally continuous and is linear inside every interval $[t_k,t_{k+1}]$ . Consequently, any linear combination of them will have the same property. Furthermore, any such function is expressible as a unique linear combination of hat functions, i.e.,

\sum_{k=0}^n c_k H_k(x)

(5.2.3)

for some choice of the coefficients $c_0,\ldots,c_n$ . No smaller set of functions can have the same properties. We summarize these facts by calling the hat functions a basis of the set of functions that are continuous and piecewise linear relative to $\mathbf{t}$ . Another point of view, familiar from abstract linear algebra, is that a basis sets up a one-to-one correspondence between the spanned function space and the more familiar space $\mathbb{R}^{n+1}$ , with each function being represented by its coefficients $c_0,\ldots,c_n$ .

An appealing characteristic of the hat function basis is that it depends only on the node locations, while the expansion coefficients in (5.2.3) depend only on the data values. This clean separation would be useful if we wanted to construct many interpolants on the same node set, and it has deeper theoretical uses as well.

Function 5.2.1 presents a simple implementation of hat functions. The inputs are a presorted vector of nodes and a value of $k$ between 0 and $n$ , which represent the indices of the endpoints. The return value is a function of $x$ that can be evaluated as needed. Note that we have not formally defined values for a hat function outside the node interval; our choice in Function 5.2.1 is to make it zero there.

Algorithm 5.2.1 (hatfun)

Julia

MATLAB

Python

Hat function

hatfun.jl

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
"""
    hatfun(t, k)

Create a piecewise linear hat function, where `t` is a
vector of n+1 interpolation nodes and `k` is an integer in 0:n
giving the index of the node where the hat function equals one.
"""
function hatfun(t, k)
    n = length(t) - 1
    return function (x)
        if k > 0 && t[k] ≤ x ≤ t[k+1]
            return (x - t[k]) / (t[k+1] - t[k])
        elseif k < n && t[k+1] ≤ x ≤ t[k+2]
            return (t[k+2] - x) / (t[k+2] - t[k+1])
        else
            return 0
        end
    end
end

Hat function

hatfun.m

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
function H = hatfun(t, k)
% HATFUN   Hat function/piecewise linear basis function.
% Input: 
%   t      interpolation nodes (vector, length n+1)
%   k      node index (integer, in 0,...,n)
% Output:
%   H      kth hat function (function)

    n = length(t) - 1;
    H = @evaluate;

    function y = evaluate(x)
        y = zeros(size(x));
        for i = 1:numel(x)
            if (k > 0) && (t(k) <= x(i)) && (x(i) <= t(k+1))
                y(i) = (x(i) - t(k)) / (t(k+1) - t(k));
            elseif (k < n) && (t(k+1) <= x(i)) && (x(i) <= t(k+2))
                y(i) = (t(k+2) - x(i)) / (t(k+2) - t(k+1));
            end
        end
    end
end

Hat function

hatfun.py

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
def hatfun(t, k):
    """
    hatfun(t, k)

    Returns a piecewise linear "hat" function,  where t is a vector of
    n+1 interpolation nodes and k is an integer in 0:n giving the index of the node
    where the hat function equals one.
    """
    n = len(t) - 1

    def evaluate(x):
        H = np.zeros(np.array(x).shape)
        for (j, xj) in enumerate(x):
            if (k > 0) and (t[k-1] <= xj) and (xj <= t[k]):
                H[j] = (xj - t[k-1]) / (t[k] - t[k-1])
            elif (k < n) and (t[k] <= xj) and (xj <= t[k+1]):
                H[j] = (t[k+1] - xj) / (t[k+1] - t[k])
        return H
    return evaluate

Example 5.2.1 (A look at hat functions)

Julia

MATLAB

Python

Example 5.2.1

Let’s define a set of four nodes (i.e., $n=3$ in our formulas).

t = [0, 0.55, 0.7, 1]

4-element Vector{Float64}:
 0.0
 0.55
 0.7
 1.0

We plot the hat functions $H_0,\ldots,H_3$ .

using Plots
plt = plot(layout=(4, 1),  legend=:top,
    xlabel=L"x",  ylims=[-0.1, 1.1],  ytick=[])
for k in 0:3
    Hₖ = FNC.hatfun(t, k)
    plot!(Hₖ, 0, 1, subplot=k + 1)
    scatter!(t, Hₖ.(t), m=3, subplot=k + 1)
    annotate!(t[k+1], 0.25, text(latexstring("H_$k"), 10), subplot=k+1)
end
plt

Example 5.2.1

Let’s define a set of four nodes (i.e., $n=3$ in our formulas).

t = [0, 0.55, 0.7, 1];

We plot the hat functions $H_0,\ldots,H_3$ .

clf
for k = 0:3
    subplot(4, 1, k+1)
    Hk = hatfun(t, k);
    fplot(Hk, [0, 1])
    hold on
    scatter(t, Hk(t))
    text(t(k+1), 0.6, sprintf("H_%d", k))
end

Unrecognized function or variable 'hatfun'.

Example 5.2.1

Let’s define a set of four nodes (i.e., $n=3$ in our formulas).

t = array([0, 0.075, 0.25, 0.55, 0.7, 1])

We plot the hat functions $H_0,\ldots,H_3$ .

x = linspace(0, 1, 300)
for k in range(6):
    plot(x, FNC.hatfun(t, k)(x))
xlabel("$x$"),  ylabel("$H_k(x)$")
title("Hat functions");

5.2.2Cardinality conditions¶

A handy property of the hat functions is that they are cardinal functions for piecewise linear interpolation, since they satisfy the cardinality conditions

H_k(t_i) = \begin{cases} 1 &\text{if $i=k$,}\\ 0 & \text{otherwise.} \end{cases}

(5.2.4)

All candidate piecewise linear (PL) functions can be expressed as a linear combination such as (5.2.3) for some coefficients $c_0,\ldots,c_n$ . But because of the cardinality conditions and the necessity for $p(x)$ to interpolate the data values in $\mathbf{y}$ , expressing the interpolant using the hat functions is trivial:

p(x) = \sum_{k=0}^n y_k H_k(x).

(5.2.5)

The resulting algorithmic simplicity is reflected in Function 5.2.2. Take note that the output of Function 5.2.2 is itself a function, meant to be called with a single argument representing a value of $x$ . Our mathematical viewpoint is that the result of an interpolation process is a function, and our codes reflect this.

Algorithm 5.2.2 (plinterp)

Julia

MATLAB

Python

Piecewise linear interpolation

plinterp.jl

1
2
3
4
5
6
7
8
9
10
11
"""
    plinterp(t, y)

Construct a piecewise linear interpolating function for data values in
`y` given at nodes in `t`.
"""
function plinterp(t, y)
    n = length(t) - 1
    H = [hatfun(t, k) for k in 0:n]
    return x -> sum(y[k+1] * H[k+1](x) for k in 0:n)
end

Piecewise linear interpolation

plinterp.m

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
function p = plinterp(t, y)
% PLINTERP   Piecewise linear interpolation.
% Input:
%   t     interpolation nodes (vector, length n+1)
%   y     interpolation values (vector, length n+1)
% Output:
%   p     piecewise linear interpolant (function)

    n = length(t) - 1;
    H = {};
    for k = 0:n
        H{k+1} = hatfun(t, k);
    end
    p = @evaluate;

    % This function evaluates p when called.
    function f = evaluate(x)
        f = 0;
        for i = 0:n
            f = f + y(i+1) * H{i+1}(x);
        end
    end
end

Piecewise linear interpolation

plinterp.py

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
def plinterp(t, y):
    """
    plinterp(t, y)

    Create a piecewise linear interpolating function for data values in y given at nodes
    in t.
    """
    n = len(t) - 1
    H = [hatfun(t, k) for k in range(n+1)]
    def evaluate(x):
        f = 0
        for k in range(n+1):
            f += y[k] * H[k](x)
        return f
    return evaluate

Example 5.2.2 (Using piecewise linear interpolation)

Julia

MATLAB

Python

Example 5.2.2

We generate a piecewise linear interpolant of $f(x)=e^{\sin 7x}$ .

f = x -> exp(sin(7x))

plot(f, 0, 1, label="function", xlabel=L"x", ylabel=L"y")

First we sample the function to create the data.

t = [0, 0.075, 0.25, 0.55, 0.7, 1]    # nodes
y = f.(t)                             # function values
scatter!(t, y, label="values at nodes")

Now we create a callable function that will evaluate the piecewise linear interpolant at any $x$ , and then plot it.

p = FNC.plinterp(t, y)
plot!(p, 0, 1, label="interpolant", title="PL interpolation")

Example 5.2.2

We generate a piecewise linear interpolant of $f(x)=e^{\sin 7x}$ .

f = @(x) exp(sin(7 * x));
clf
fplot(f, [0, 1], displayname="function")
xlabel("x");  ylabel(("y"));

First we sample the function to create the data.

t = [0, 0.075, 0.25, 0.55, 0.7, 1];    % nodes
y = f(t);                              % function values

Now we create a callable function that will evaluate the piecewise linear interpolant at any $x$ , and then plot it.

p = plinterp(t, y);
hold on
fplot(p, [0, 1], displayname="interpolant")
scatter(t, y, displayname="values at nodes")
title("PL interpolation")
legend();

Unrecognized function or variable 'plinterp'.

Example 5.2.2

We generate a piecewise linear interpolant of $f(x)=e^{\sin 7x}$ .

f = lambda x: exp(sin(7 * x))
x = linspace(0, 1, 400)
fig, ax = subplots()
plot(x, f(x), label="function")
xlabel("$x$")
ylabel("$f(x)$");

First we sample the function to create the data.

t = array([0, 0.075, 0.25, 0.55, 0.7, 1])  # nodes
y = f(t)  # function values

ax.plot(t, y, "o", label="nodes")
ax.legend()
fig

Now we create a callable function that will evaluate the piecewise linear interpolant at any $x$ , and then plot it.

p = FNC.plinterp(t, y)
ax.plot(x, p(x), label="interpolant")
ax.legend()
fig

5.2.3Conditioning and convergence¶

The condition number bounds from Theorem 5.1.1 are very simple for piecewise linear interpolation because the interpolant of the data $\mathbf{e}_k$ is just the hat function $H_k$ . Hence, $1\le \kappa \le n+1$ . However, there is an even simpler result.

Proof

By linearity,

\mathcal{I}(\mathbf{y}+\mathbf{z}) - \mathcal{I}(\mathbf{y}) = \mathcal{I}(\mathbf{z}) = \sum_{k=0}^n z_k H_k(x).

(5.2.7)

Call this piecewise linear function $p(x)$ . Consider a maximum element of $\mathbf{z}$ , i.e., choose $i$ such that $|z_i|=\|\mathbf{z}\|_\infty$ . Then $|p(t_i)|=\|\mathbf{z}\|_\infty$ . Hence $\|p\|_\infty\ge \|\mathbf{z}\|_\infty$ . Now consider

|p(x)| = \left|\sum_{k=0}^n z_k H_k(x)\right| \le \sum_{k=0}^n |z_k| H_k(x) \le \|\mathbf{z}\|_\infty \sum_{k=0}^n H_k(x) = \|\mathbf{z}\|_\infty.

(5.2.8)

You are asked to prove the final step above in Exercise 5.2.4. We conclude that $\|p\|_\infty\le \|\mathbf{z}\|_\infty$ , so that $\|p\|_\infty = \|\mathbf{z}\|_\infty$ , which completes the proof.

Now suppose that $f$ is a “nice” function on an interval $[a,b]$ containing all the nodes. We can sample values of $f$ to get data, i.e., $y_k=f(t_k)$ for all $k$ , then perform piecewise linear interpolation of the data to get a different function, the interpolant $p$ . How close is $p$ to the original $f$ ?

To make a simple statement, we will consider only the case of equally spaced nodes covering the interval. It turns out that piecewise linear interpolation converges at second order in the spacing of the nodes.

For an outline of a proof, see Exercise 5.2.5.

We normally don’t have access to $f''$ , so the importance of Theorem 5.2.2 is that the error in the interpolant is $O(h^2)$ as $h\to 0$ .

Thus, Theorem 5.2.2 states that piecewise linear interpolation is second-order accurate. For instance, if we increase the number of equally spaced nodes by a factor of 10, the piecewise linear interpolant becomes about 100 times more accurate. Note also that if $y \approx C h^m$ , then

\log y \approx m (\log h) + \log C.

(5.2.10)

Example 5.2.3 (Convergence of piecewise linear interpolation)

Julia

MATLAB

Python

Example 5.2.3

We measure the convergence rate for piecewise linear interpolation of $e^{\sin 7x}$ over $x \in [0,1]$ .

f = x -> exp(sin(7x))
x = range(0, 1, 10001)  # sample the difference at many points
n = @. round(Int, 10^(1:0.25:3.5))
maxerr = zeros(length(n))
for (k, n) in enumerate(n)
    t = (0:n) / n    # interpolation nodes
    p = FNC.plinterp(t, f.(t))
    err = @. f(x) - p(x)
    maxerr[k] = norm(err, Inf)
end

data = (n=n[1:4:end], err=maxerr[1:4:end])
@pt :header=["n", "max-norm error"] data

As predicted, a factor of 10 in $n$ produces a factor of 100 in the error. In a convergence plot, it is traditional to have $h$ decrease from left to right, so we expect a straight line of slope -2 on a log-log plot.

h = @. 1 / n
order2 = @. 10 * (h / h[1])^2

plot(h, maxerr, m=:o, label="error", xflip=true)
plot!(h, order2;
    l=:dash,  label=L"O(h^2)",
    xaxis=(:log10, L"h"),  yaxis=(:log10, L"|| f-p\, ||_\infty"),
    title="Convergence of PL interpolation")

Example 5.2.3

We measure the convergence rate for piecewise linear interpolation of $e^{\sin 7x}$ over $x \in [0,1]$ .

f = @(x) exp(sin(7 * x));
x = linspace(0, 1, 10001)';    % sample the difference at many points
n = round(10.^(1:0.25:3.5))';
maxerr = zeros(size(n));
for i = 1:length(n)
    t = (0:n(i)) / n(i);       % interpolation nodes
    p = plinterp(t, f(t));
    maxerr(i) = norm(f(x) - p(x), Inf);
end
table(n(1:4:end), maxerr(1:4:end), variableNames=["n", "inf-norm error"])

Unrecognized function or variable 'plinterp'.

As predicted, a factor of 10 in $n$ produces a factor of 100 reduction in the error. In a convergence plot, it is traditional to have $h$ decrease from left to right, so we expect a straight line of slope -2 on a log-log plot.

clf
loglog(n, maxerr, "-o", displayname="error")
order2 = 0.5 * maxerr(end) * (n / n(end)) .^ (-2);
hold on
loglog(n, order2, "k--", displayname="O(n^{-2})")
xlabel("n");  ylabel("|| f-p ||_{\infty}")
title("Convergence of PL interpolation")
legend();

Example 5.2.3

We measure the convergence rate for piecewise linear interpolation of $e^{\sin 7x}$ over $x \in [0,1]$ .

f = lambda x: exp(sin(7 * x))
x = linspace(0, 1, 10000)  # sample the difference at many points
N = 2 ** arange(3, 11)
err = zeros(N.size)
for i, n in enumerate(N):
    t = linspace(0, 1, n + 1)  # interpolation nodes
    p = FNC.plinterp(t, f(t))
    err[i] = max(abs(f(x) - p(x)))
print(err)

[2.16029984e-01 6.38173511e-02 1.60381329e-02 4.05882168e-03
 1.01556687e-03 2.54022468e-04 6.35007579e-05 1.58778800e-05]

order2 = 0.1 * (N / N[0]) ** (-2)
loglog(N, err, "-o", label="observed error")
loglog(N, order2, "--", label="2nd order")
xlabel("$n$")
ylabel("$\|f-p\|_\infty$")
legend();

<>:5: SyntaxWarning: invalid escape sequence '\|'
<>:5: SyntaxWarning: invalid escape sequence '\|'
/var/folders/0m/fy_f4_rx5xv68sdkrpcm26r00000gq/T/ipykernel_23602/3346807096.py:5: SyntaxWarning: invalid escape sequence '\|'
  ylabel("$\|f-p\|_\infty$")

5.2.4Exercises¶

Exercise 5.2.1

⌨ For each given function and interval, perform piecewise linear interpolation using Function 5.2.2 for $n+1$ equispaced nodes with $n=10,20,40,80,160,320$ . For each $n$ , estimate the error

E(n) = \| f-p \|_\infty = \max_x | f(x) - p(x) |

(5.2.11)

by evaluating the function and interpolant at 1600 points in the interval. Make a log-log plot of $E$ as a function of $n$ and add the line $E=Cn^{-2}$ for a constant $C$ of your choosing.

(a) $\cos(\pi x^2)$ on $[0,4]$

(b) $\log(x)$ on $[1,20]$

(c) $\sin\left(\frac{1}{x}\right)$ on $\left[\frac{1}{2},7\right]$

Exercise 5.2.5

✍ Here we consider a proof of Theorem 5.2.2 using the mean value theorems from elementary calculus: If $f$ is continuously differentiable in $(a,b)$ , then there exist points $s$ and $t$ in $(a,b)$ such that

\int_a^b f(z) \, dz = (b-a)f(s) \qquad \text{and} \qquad f'(t) = \frac{f(b)-f(a)}{b-a}.

(5.2.13)

For the following, suppose $x \in (t_k,t_{k+1})$ .

(a) Show that for some $s \in (t_k,t_{k+1})$ ,

f(x) = y_k + (x-t_k)f'(s).

(5.2.14)

(b) Show that for some other values $u$ and $v$ in $(t_k,t_{k+1})$ ,

f'(s) - \frac{y_{k+1}-y_k}{t_{k+1}-t_k} = (s-u) f''(v).

(5.2.15)

(c) Use (5.2.1) to finish the proof of the theorem.

Preface

The interpolation problem

Preface

Cubic splines