The barycentric formula Tobin A. Driscoll Richard J. Braun
The Lagrange formula (9.1.4) x x x β k \ell_k β k β x x x n n n O ( n 2 ) O(n^2) O ( n 2 ) x x x
9.2.1 Derivation ΒΆ We again will use the error indicator function Ξ¦ \Phi Ξ¦ Definition 9.1.2
Ξ¦ ( x ) = β j = 0 n ( x β t j ) , \Phi(x) = \prod_{j=0}^n (x-t_j), Ξ¦ ( x ) = j = 0 β n β ( x β t j β ) , as well as a set of values derived from the nodes.
The barycentric weights for nodes x 0 , β¦ , x n x_0,\dots,x_n x 0 β , β¦ , x n β
w k = 1 β j = 0 j β k n ( t k β t j ) = 1 Ξ¦ β² ( t k ) , k = 0 , β¦ , n . w_k = \frac{1}{\displaystyle \prod_{\substack{j=0\\j\neq k}}^n (t_k - t_j)} = \frac{1}{\Phi'(t_k)}, \qquad
k = 0,\ldots,n. w k β = j = 0 j ξ = k β β n β ( t k β β t j β ) 1 β = Ξ¦ β² ( t k β ) 1 β , k = 0 , β¦ , n . The following formula is the key to efficient and stable evaluation of a polynomial interpolant.
The Lagrange cardinal polynomial (9.1.3)
β k ( x ) = Ξ¦ ( x ) w k x β t k , \ell_k(x) = \Phi(x) \frac{w_k}{x-t_k}, β k β ( x ) = Ξ¦ ( x ) x β t k β w k β β , and thus the interpolating polynomial in (9.1.4)
p ( x ) = Ξ¦ ( x ) β k = 0 n w k x β t k y k . p(x) = \Phi(x) \sum_{k=0}^n \frac{w_k}{x-t_k} y_k. p ( x ) = Ξ¦ ( x ) k = 0 β n β x β t k β w k β β y k β . Obviously, the constant function p ( x ) β‘ 1 p(x)\equiv 1 p ( x ) β‘ 1 TheoremΒ 9.1.1 y k = 1 y_k=1 y k β = 1 k k k (9.2.5)
1 = Ξ¦ ( x ) β k = 0 n w k x β t k . 1 = \Phi(x) \sum_{k=0}^n \frac{w_k}{x-t_k}. 1 = Ξ¦ ( x ) k = 0 β n β x β t k β w k β β . This is solved for Ξ¦ ( x ) \Phi(x) Ξ¦ ( x ) (9.2.5) (9.2.3)
Equation (9.2.3) x x x lim β‘ x β t k p ( x ) = y k \lim_{x\to t_k} p(x) = y_k lim x β t k β β p ( x ) = y k β ExerciseΒ 9.2.3
Let us write out the barycentric formula for the interpolating polynomial for the quadratic case (n = 2 n=2 n = 2 Example 9.1.2 (9.2.2)
w 0 = 1 ( t 0 β t 1 ) ( t 0 β t 2 ) = 1 ( 0 β Ο 6 ) ( 0 β Ο 3 ) = 18 Ο 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}, w 0 β = ( t 0 β β t 1 β ) ( t 0 β β t 2 β ) 1 β = ( 0 β 6 Ο β ) ( 0 β 3 Ο β ) 1 β = Ο 2 18 β , and, similarly, w 1 = β 36 / Ο 2 w_1 = -36/\pi^2 w 1 β = β 36/ Ο 2 w 2 = 18 / Ο 2 w_2=18/\pi^2 w 2 β = 18/ Ο 2
Note that in (9.2.3) w 0 = w 2 = 1 w_0=w_2=1 w 0 β = w 2 β = 1 w 1 = β 2 w_1=-2 w 1 β = β 2
p ( x ) = w 0 x β t 0 y 0 + w 1 x β t 1 y 1 + w 2 x β t 2 y 2 w 0 x β t 0 + w 1 x β t 1 + w 2 x β t 2 = ( 1 x ) 0 β ( 2 x β Ο / 6 ) 1 3 + ( 1 x β Ο / 3 ) 3 1 x β 2 x β Ο / 6 + 1 x β Ο / 3 . \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*} p ( x ) β = x β t 0 β w 0 β β + x β t 1 β w 1 β β + x β t 2 β w 2 β β x β t 0 β w 0 β β y 0 β + x β t 1 β w 1 β β y 1 β + x β t 2 β w 2 β β y 2 β β = x 1 β β x β Ο /6 2 β + x β Ο /3 1 β ( x 1 β ) 0 β ( x β Ο /6 2 β ) 3 β 1 β + ( x β Ο /3 1 β ) 3 β β . β Further algebraic manipulation could return this expression to the classical Lagrange form derived in Example 9.1.2
9.2.2 Implementation ΒΆ For certain important node distributions, simple formulas for the weights w k w_k w k β w k w_k w k β w k β 1 w_k^{-1} w k β 1 β
We begin with the singleton node set { t 0 } \{t_0\} { t 0 β } w 0 = 1 w_0=1 w 0 β = 1 Ο k , m β 1 \omega_{k,m-1} Ο k , m β 1 β k < m k< m k < m k k k { t 0 , β¦ , t m β 1 } \{t_0,\ldots,t_{m-1}\} { t 0 β , β¦ , t m β 1 β }
Ο k , m = β j = 0 j β k m ( t k β t j ) = Ο k , m β 1 β
( t k β t m ) , k = 0 , 1 , β¦ , m β 1. \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. Ο k , m β = j = 0 j ξ = k β β m β ( t k β β t j β ) = Ο k , m β 1 β β
( t k β β t m β ) , k = 0 , 1 , β¦ , m β 1. A direct application of (9.2.2) Ο m , m \omega_{m,m} Ο m , m β m = 1 , β¦ , n m=1,\ldots,n m = 1 , β¦ , n w k = Ο k , n β 1 w_k=\omega_{k,n}^{-1} w k β = Ο k , n β 1 β
In Function 9.2.1
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
"""
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)
Ξ = x .- t
if any(iszero.(Ξ)) # we're at a node exactly
# return the node's data value
idx = findfirst(iszero.(Ξ))
f = y[idx]
else
terms = w ./ Ξ
f = sum(y .* terms) / sum(terms)
end
end
return p
end
As noted in Example 9.2.1 (9.2.3)
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 x x
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
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(:)))' % where divided by zero
[~, idx] = min(abs(x(j) - t)); % node closest to x(j)
f(j) = y(idx); % value at node
end
end % evaluate function
end % main function1
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)
As noted in Example 9.2.1 (9.2.3)
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 x x
Computing all n + 1 n+1 n + 1 Function 9.2.1 O ( n 2 ) O(n^2) O ( n 2 ) p ( x ) p(x) p ( x ) x x x O ( n ) O(n) O ( n )
We show the barycentric formula in action for values from the function sin β‘ ( e 2 x ) \sin(e^{2x}) sin ( e 2 x ) [ 0 , 1 ] [0,1] [ 0 , 1 ] n = 3 n=3 n = 3 n = 6 n=6 n = 6
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 n=6 n = 6
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")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 n=3 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 n=6 n = 6
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")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")
figThe curves must intersect at the interpolation nodes. For n = 6 n=6 n = 6
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.3 Stability ΒΆ You might suspect that as the evaluation point x x x t k t_k t k β 1 / ( x β t k ) 1/(x-t_k) 1/ ( x β t k β ) 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.4 Exercises ΒΆ β (a) Find the barycentric weights for the nodes t 0 = 0 t_0=0 t 0 β = 0 t 1 = 1 t_1=1 t 1 β = 1 t 2 = 3 t_2=3 t 2 β = 3
(b) Compute the interpolant at x = 2 x=2 x = 2 y 0 = β 2 y_0=-2 y 0 β = β 2 y 1 = 2 y_1=2 y 1 β = 2 y 2 = 1 y_2=1 y 2 β = 1
β Show using LβHΓ΄pitalβs rule on (9.2.3) p ( t i ) = y i p(t_i)=y_i p ( t i β ) = y i β i = 0 , β¦ , n i=0,\ldots,n i = 0 , β¦ , n
β¨ In each case, use Function 9.2.1 n + 1 n+1 n + 1 n n n
(a) f ( x ) = ln β‘ ( x ) , n = 2 , 3 , 4 , x β [ 1 , 10 ] f(x) = \ln (x), \quad n = 2,3,4, \quad x\in [1,10] f ( x ) = ln ( x ) , n = 2 , 3 , 4 , x β [ 1 , 10 ]
(b) f ( x ) = tanh β‘ ( x ) , n = 2 , 3 , 4 , x β [ β 3 , 2 ] f(x) = \tanh (x), \quad n = 2,3,4, \quad x \in [-3,2] f ( x ) = tanh ( x ) , n = 2 , 3 , 4 , x β [ β 3 , 2 ]
(c) f ( x ) = cosh β‘ ( x ) , n = 2 , 3 , 4 , x β [ β 1 , 3 ] f(x) = \cosh (x), \quad n = 2,3,4, \quad x \in [-1,3] f ( x ) = cosh ( x ) , n = 2 , 3 , 4 , x β [ β 1 , 3 ]
(d) f ( x ) = β£ x β£ , n = 3 , 5 , 7 , x β [ β 2 , 1 ] f(x) = |x|, \quad n = 3,5,7, \quad x \in [-2,1] f ( x ) = β£ x β£ , n = 3 , 5 , 7 , x β [ β 2 , 1 ]
β¨ Using code from Function 9.2.1 n + 1 n+1 n + 1 [ β 1 , 1 ] [-1,1] [ β 1 , 1 ] n = 30 n=30 n = 30 n = 60 n=60 n = 60 n = 90 n=90 n = 90 β£ w i β£ |w_i| β£ w i β β£ t i t_i t i β Stability of polynomial interpolation .)