5392320Interactive illustration of Floater-Hormann barycentric rational interpolationJeremy Magland 1% ---------------------------------------------------------------------------
2% The polynomial interpolant, in barycentric form: the d = n member
3%
4% Equation (2), first written down by Taylor and by Dupuy:
5%
6% w_k = prod_{j ~= k} 1 / (x_k - x_j).
7%
8% These are what equation (18) gives when d = n, since then there is a single
9% local polynomial and it is the interpolating polynomial p_n itself. The
10% paper's remark that the weights of (2) "prevent poles" is visible on the
11% Poles tab in a degenerate way: the denominator of r reduces to the constant 1,
12% because the Lagrange basis functions sum to 1, so there is nothing to vanish
13% and no roots at all.
14%
15% Having no poles is not the same as approximating well. Leave the function at
16% Runge's 1/(1+x^2), the nodes uniform, and push n up: this is the divergence
17% the paper opens with. Switch the nodes to Chebyshev and it behaves.
18%
19% The d slider is ignored by this script.
20% ---------------------------------------------------------------------------
22function w = bary_weights(x, d)
23n1 = numel(x);
24w = zeros(1, n1);
25for k = 1:n1
26 p = 1;
27 for j = 1:n1
28 if j ~= k
29 p = p / (x(k) - x(j));
30 end
31 end
32 w(k) = p;
33end
34end
36function r = bary_eval(x, y, w, t)
37sz = size(t);
38D = t(:) - x(:).';
39Q = w(:).' ./ D;
40r = (Q * y(:)) ./ sum(Q, 2);
41hit = find(any(D == 0, 2));
42for m = 1:numel(hit)
43 k = find(D(hit(m), :) == 0, 1);
44 r(hit(m)) = y(k);
45end
46r = reshape(r, sz);
47end