1% ---------------------------------------------------------------------------
2% Berrut's interpolant: the d = 0 member of the family
3% J.-P. Berrut, Comput. Math. Appl. 15 (1988) 1-16
4%
5% Equation (3) of the paper. The weights are simply the alternating signs,
6%
7% w_k = (-1)^k,
8%
9% which is what equation (18) reduces to when d = 0 (up to a common positive
10% factor, which does not change r). Berrut showed this has no real poles;
11% Floater and Hormann's Theorem 3 gives it the rate O(h), but only under a
12% bound on the local mesh ratio beta. Set the nodes to "paired" and watch what
13% happens: the interpolant develops kinks and the error stops falling, while
14% the d >= 1 members of the family are untroubled.
15%
16% The d slider is ignored by this script.
17% ---------------------------------------------------------------------------
19function w = bary_weights(x, d)
20n = numel(x) - 1;
21w = (-1).^(0:n);
22end
24function r = bary_eval(x, y, w, t)
25sz = size(t);
26D = t(:) - x(:).';
27Q = w(:).' ./ D;
28r = (Q * y(:)) ./ sum(Q, 2);
29hit = find(any(D == 0, 2));
30for m = 1:numel(hit)
31 k = find(D(hit(m), :) == 0, 1);
32 r(hit(m)) = y(k);
33end
34r = reshape(r, sz);
35end
37function [P, L] = local_blend(x, y, d, t)
38% With d = 0 the "local polynomials" are the constants p_i = f_i, and the
39% blending functions are mu_i(t) = prod_{j<i} (t - x_j) * prod_{k>i} (x_k - t)
40% normalised to sum to 1.
41n = numel(x) - 1;
42t = reshape(t, 1, []);
43P = repmat(y(:), 1, numel(t));
44Mu = zeros(n + 1, numel(t));
45for i = 0:n
46 pr = ones(1, numel(t));
47 for j = 0:(i - 1)
48 pr = pr .* (t - x(j + 1));
49 end
50 for k = (i + 1):n
51 pr = pr .* (x(k + 1) - t);
52 end
53 Mu(i + 1, :) = pr;
54end
55L = Mu ./ sum(Mu, 1);
56end