from numpy import *
from scipy import linalg
from scipy.linalg import norm
from matplotlib.pyplot import *
from prettytable import PrettyTable
from timeit import default_timer as timer
import sys
sys.path.append('fncbook/')
import fncbook as FNC

# This (optional) block is for improving the display of plots.
# from IPython.display import set_matplotlib_formats
# set_matplotlib_formats("svg","pdf")
# %config InlineBackend.figure_format = 'svg'
rcParams["figure.figsize"] = [7, 4]
rcParams["lines.linewidth"] = 2
rcParams["lines.markersize"] = 4
rcParams['animation.html'] = "jshtml"  # or try "html5"

9.2.The barycentric formula¶

The Lagrange formula (9.1.4) is useful theoretically but not ideal for computation. For each new value of $x$ , all 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 $\Phi$ 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.

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 9.2.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}}\\[2ex] & =\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*}

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.8)

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)

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)

About the code

As noted in Example 9.2.1, a common scaling factor in the weights does not affect the barycentric formula (9.2.3). In lines 9--10 this fact is used to rescale the nodes in order to avoid eventual tiny or enormous numbers that could go outside the bounds of double precision.

The return value is a function that evaluates the polynomial interpolant. Within this function, isinf is used to detect either Inf or -Inf, which occurs when $x$ exactly equals one of the nodes. In this event, the corresponding data value is returned.

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$ .

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.

The barycentric formula

9.2.The barycentric formula¶

9.2.1Derivation¶

9.2.2Implementation¶

9.2.3Stability¶

9.2.4Exercises¶