The barycentric formula - Fundamentals of Numerical Computation

The Lagrange formula (9.1.4) is useful theoretically but not ideal for computation. For each new value of $x$ , all of the cardinal functions $\ell_k$ must be evaluated at $x$ , which requires a product of $n$ terms. Thus the total work is $O(n^2)$ for every value of $x$ . Moreover, the formula is numerically unstable. An alternative version of the formula improves on both issues.

9.2.1Derivation¶

We again will use the error indicator function Φ from Definition 9.1.2,

\Phi(x) = \prod_{j=0}^n (x-t_j),

(9.2.1)

as well as a set of values derived from the nodes.

The following formula is the key to efficient and stable evaluation of a polynomial interpolant.

Proof 9.2.1

The Lagrange cardinal polynomial (9.1.3) can be written as

\ell_k(x) = \Phi(x) \frac{w_k}{x-t_k},

(9.2.4)

and thus the interpolating polynomial in (9.1.4) is

p(x) = \Phi(x) \sum_{k=0}^n \frac{w_k}{x-t_k} y_k.

(9.2.5)

Obviously, the constant function $p(x)\equiv 1$ is its own polynomial interpolant on any set of nodes. The uniqueness of the interpolating polynomial, as proved in Theorem 9.1.1, allows us to plug $y_k=1$ for all $k$ into (9.2.5) to obtain

1 = \Phi(x) \sum_{k=0}^n \frac{w_k}{x-t_k}.

(9.2.6)

This is solved for $\Phi(x)$ and put back into (9.2.5) to get (9.2.3).

Equation (9.2.3) is certainly an odd-looking way to write a polynomial! Indeed, it is technically undefined when $x$ equals one of the nodes, but in fact, $\lim_{x\to t_k} p(x) = y_k$ , so a continuous extension to the nodes is justified. (See Exercise 3.)

Example 9.2.1

Let us write out the barycentric formula for the interpolating polynomial for the quadratic case ( $n=2$ ) for Example 9.1.2. The weights are computed from (9.2.2):

w_0 = \frac{1}{(t_0-t_1)(t_0-t_2)} = \frac{1}{\left(0-\frac{\pi}{6}\right) \left(0-\frac{\pi}{3}\right)} = \frac{18}{\pi^2},

(9.2.7)

and similarly, $w_1 = -36/\pi^2$ and $w_2=18/\pi^2$ .

Note that in (9.2.3), any common factor in the weights cancels out without affecting the results. Hence it’s a lot easier to use $w_0=w_2=1$ and $w_1=-2$ . Then

\begin{align*} p(x) & = \frac{\rule[-1.2em]{0pt}{1em} \dfrac{w_0}{x-t_0} y_0 + \dfrac{w_1}{x-t_1} y_1 + \dfrac{w_2}{x-t_2} y_2 }{ \rule{0pt}{1.5em} \dfrac{w_0}{x-t_0} + \dfrac{w_1}{x-t_1} + \dfrac{w_2}{x-t_2}}\\[1.5ex] & =\frac{ \rule[-1.2em]{0pt}{1em}\left( \dfrac{1}{x} \right) 0 - \left( \dfrac{2}{x-\pi/6} \right) \dfrac{1}{\sqrt{3}} + \left( \dfrac{1}{x-\pi/3} \right) \sqrt{3} }{ \rule{0pt}{1.6em} \dfrac{1}{x} - \dfrac{2}{x-\pi/6} + \dfrac{1}{x-\pi/3} }. \end{align*}

(9.2.8)

Further algebraic manipulation could return this expression to the classical Lagrange form derived in Example 9.1.2.

9.2.2Implementation¶

For certain important node distributions, simple formulas for the weights $w_k$ are known. Otherwise, the first task of an implementation is to compute the weights $w_k$ , or more conveniently, $w_k^{-1}$ .

We begin with the singleton node set $\{t_0\}$ , for which one gets the single weight $w_0=1$ . The idea is to grow this singleton into the set of all nodes through a recursive formula. Define $\omega_{k,m-1}$ (for $k< m$ ) as the inverse of the weight for node $k$ using the set $\{t_0,\ldots,t_{m-1}\}$ . Then

\omega_{k,m} = \displaystyle \prod_{\substack{j=0\\j\neq k}}^{m} (t_k - t_j) = \omega_{k,m-1} \cdot (t_k-t_{m}), \qquad k=0,1,\ldots,m-1.

(9.2.9)

A direct application of (9.2.2) can be used to find $\omega_{m,m}$ . This process is iterated over $m=1,\ldots,n$ to find $w_k=\omega_{k,n}^{-1}$ .

In Function 9.2.1 we give an implementation of the barycentric formula for polynomial interpolation.

Algorithm 9.2.1 (polyinterp)

Julia

MATLAB

Python

Barycentric polynomial interpolation

polyinterp.jl

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
"""
    polyinterp(t, y)

Construct a callable polynomial interpolant through the points in
vectors `t`, `y` using the barycentric interpolation formula.
"""
function polyinterp(t, y)
    n = length(t) - 1
    C = (t[n+1] - t[1]) / 4           # scaling factor to ensure stability
    tc = t / C

    # Adding one node at a time, compute inverses of the weights.
    ω = ones(n + 1)
    for m in 0:n-1
        d = tc[1:m+1] .- tc[m+2]    # vector of node differences
        @. ω[1:m+1] *= d            # update previous
        ω[m+2] = prod(-d)         # compute the new one
    end
    w = 1 ./ ω                      # go from inverses to weights

    # This function evaluates the interpolant at given x.
    p = function (x)
        terms = @. w / (x - t)
        if any(isinf.(terms))     # there was division by zero
            # return the node's data value
            idx = findfirst(x .== t)
            f = y[idx]
        else
            f = sum(y .* terms) / sum(terms)
        end
    end
    return p
end

Barycentric polynomial interpolation

polyinterp.m

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
function p = polyinterp(t,y)
% POLYINTERP Polynomial interpolation by the barycentric formula.
% Input:
%   t   interpolation nodes (vector, length n+1)
%   y   interpolation values (vector, length n+1)
% Output:
%   p   polynomial interpolant (function)

t = t(:);                    % column vector
n = length(t)-1;
C = (t(end)-t(1)) / 4;       % scaling factor to ensure stability
tc = t/C;

% Adding one node at a time, compute inverses of the weights.
omega = ones(n+1,1);
for m = 1:n
    d = (tc(1:m) - tc(m+1));    % vector of node differences
    omega(1:m) = omega(1:m).*d; % update previous 
    omega(m+1) = prod( -d );    % compute the new one
end
w = 1./omega;                   % go from inverses to weights

p = @evaluate;

    function f = evaluate(x)
        % % Compute interpolant, one value of x at a time.
        f = zeros(size(x));
        for j = 1:numel(x)
            terms = w ./ (x(j) - t );
            f(j) = sum(y.*terms) / sum(terms);
        end
        
        % Apply L'Hopital's Rule exactly.
        for j = find( isnan(f(:)) )'      % divided by zero here
            [~,idx] = min( abs(x(j)-t) );   % node closest to x(j)
            f(j) = y(idx);                  % value at node
        end        
    end

end

Barycentric polynomial interpolation

polyinterp.py

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
def polyinterp(t, y):
    """
    polyinterp(t, y)

    Return a callable polynomial interpolant through the points in vectors t, y. Uses
    the barycentric interpolation formula.
    """
    n = len(t) - 1
    C = (t[-1] - t[0]) / 4  # scaling factor to ensure stability
    tc = t / C

    # Adding one node at a time, compute inverses of the weights.
    omega = np.ones(n + 1)
    for m in range(n):
        d = tc[: m + 1] - tc[m + 1]  # vector of node differences
        omega[: m + 1] = omega[: m + 1] * d  # update previous
        omega[m + 1] = np.prod(-d)  # compute the new one
    w = 1 / omega  # go from inverses to weights

    def p(x):
        # Compute interpolant.
        z = np.where(x == t)[0]
        if len(z) > 0:  # avoid dividing by zero
            # Apply L'Hopital's Rule exactly.
            f = y[z[0]]
        else:
            terms = w / (x - t)
            f = np.sum(y * terms) / np.sum(terms)
        return f

    return np.vectorize(p)

Computing all $n+1$ weights in Function 9.2.1 takes $O(n^2)$ operations. Fortunately, the weights depend only on the nodes, not the data, and once they are known, computing $p(x)$ at a particular value of $x$ takes just $O(n)$ operations.

Example 9.2.2 (Barycentric interpolation)

We show the barycentric formula in action for values from the function $\sin(e^{2x})$ at equally spaced nodes in $[0,1]$ with $n=3$ and $n=6$ .

Julia

MATLAB

Python

Example 9.2.2

using Plots
f(x) = sin(exp(2x))
plot(f, 0, 1, label="function", legend=:bottomleft)

t = (0:3) / 3
y = f.(t)
scatter!(t, y, color=:black, label="nodes")

p = FNC.polyinterp(t, y)
plot!(p, 0, 1, label="interpolant", title="Interpolation on 4 nodes")

The curves must intersect at the interpolation nodes. For $n=6$ the interpolant is noticeably better.

plot(f, 0, 1, label="function", legend=:bottomleft)
t = (0:6) / 6
y = f.(t)
p = FNC.polyinterp(t, y)
scatter!(t, y, color=:black, label="nodes")
plot!(p, 0, 1, label="interpolant", title="Interpolation on 7 nodes")

Example 9.2.2

f = @(x) sin( exp(2 * x) );
clf,  fplot(f, [0, 1], displayname="function")
xlabel('x'),  ylabel('f(x)')   
legend(location="southwest");

We start with 4 equally spaced nodes ( $n=3$ ).

t = linspace(0, 1, 4)'; 
y = f(t);
p = polyinterp(t, y);
hold on,  fplot(p, [0, 1], displayname="interpolant on 4 nodes")
scatter(t, y, 'k', displayname="nodes")

The curves always intersect at the interpolation nodes. For $n=6$ , the interpolant is noticeably better.

cla,  fplot(f, [0, 1], displayname="function")
t = linspace(0, 1, 7)'; 
y = f(t);
p = polyinterp(t, y);
hold on,  fplot(p, [0, 1], displayname="interpolant on 7 nodes")
scatter(t, y, 'k', displayname="nodes")

Example 9.2.2

f = lambda x: sin(exp(2 * x))
x = linspace(0, 1, 500)
fig, ax = subplots()
ax.plot(x, f(x), label="function")

t = linspace(0, 1, 4)
y = f(t)
p = FNC.polyinterp(t, y)

ax.plot(x, p(x), label="interpolant")
ax.plot(t, y, "ko", label="nodes")
ax.legend()
ax.set_title("Interpolation on 4 nodes")
fig

The curves must intersect at the interpolation nodes. For $n=6$ the interpolant is noticeably better.

plot(x, f(x), label="function")
t = linspace(0, 1, 7)
y = f(t)
p = FNC.polyinterp(t, y)
plot(x, p(x), label="interpolant")
plot(t, y, "ko", label="nodes")
legend(),  title("Interpolation on 7 nodes");

9.2.3Stability¶

You might suspect that as the evaluation point $x$ approaches a node $t_k$ , subtractive cancellation error will creep into the barycentric formula because of the term $1/(x-t_k)$ . While such errors do occur, they turn out not to cause trouble, because the same cancellation happens in the numerator and denominator. In fact, the stability of the barycentric formula has been proved, though we do not give the details.

9.2.4Exercises¶

✍ (a) Find the barycentric weights for the nodes $t_0=0$ , $t_1=1$ , $t_2=3$ .
(b) Compute the interpolant at $x=2$ for the nodes in part (a) and the data $y_0=-2$ , $y_1=2$ , $y_2=1$ .
✍ For each case of Exercise 9.1.1, write out the barycentric form of the interpolating polynomial.

✍ Show using L’Hôpital’s rule on (9.2.3) that $p(t_i)=y_i$ for all $i=0,\ldots,n$ .
⌨ In each case, use Function 9.2.1 to interpolate the given function using $n+1$ evenly spaced nodes in the given interval. Plot each interpolant together with the exact function.
(a) $f(x) = \ln (x), \quad n = 2,3,4, \quad x\in [1,10]$
(b) $f(x) = \tanh (x), \quad n = 2,3,4, \quad x \in [0-3,2]$
(c) $f(x) = \cosh (x), \quad n = 2,3,4, \quad x \in [-1,3]$
(d) $f(x) = |x|, \quad n = 3,5,7, \quad x \in [-2,1]$
⌨ Using code from Function 9.2.1, compute the barycentric weights numerically using $n+1$ equally spaced nodes in $[-1,1]$ for $n=30$ , $n=60$ , and $n=90$ . On a single graph, plot $|w_i|$ as a function of $t_i$ on a log-linear scale. (The resulting graphs are an indication of the trouble with equally spaced nodes that is explored in Stability of polynomial interpolation.)
✍ Derive this fact stated implicitly in (9.2.2):
$\Phi'(x_k) = \prod_{\substack{j=0\\j\neq k}}^n (t_k - t_j).$
(9.2.10)
✍ Use (9.2.4) to show that if $j\neq k$ ,
$\ell_k'(x_j) = \frac{w_k}{w_j(x_j-x_k)}.$
(9.2.11)

Preface

Polynomial interpolation

Preface

Stability of polynomial interpolation