5392320Interactive illustration of Floater-Hormann barycentric rational interpolationJeremy Magland 1% ---------------------------------------------------------------------------
2% The integer weights of Section 4, for equally spaced nodes
3%
4% When the nodes are uniform with spacing h, equation (18) collapses to
5%
6% w_k = (-1)^(k-d) / h^d * sum_{i in J_k} 1 / ((k-i)! (i+d-k)!),
7%
8% and since a common positive factor does not change r we may multiply by
9% d! h^d and read off integers:
10%
11% w_k = (-1)^(k-d) * sum_{i in J_k} binomial(d, k-i).
12%
13% Writing delta_k = |w_k|, the first few rows are the ones tabulated in the
14% paper:
15%
16% d = 0: 1, 1, ..., 1, 1
17% d = 1: 1, 2, 2, ..., 2, 2, 1
18% d = 2: 1, 3, 4, ..., 4, 3, 1
19% d = 3: 1, 4, 7, 8, 8, ..., 8, 8, 7, 4, 1
20% d = 4: 1, 5, 11, 15, 16, 16, ..., 16, 16, 15, 11, 5, 1
21%
22% Almost every weight is the same; the only difference is at the two ends.
23% Yet that small change is what raises the approximation order from O(h) to
24% O(h^(d+1)). The "Blending & weights" tab checks these against the general
25% formula, so switching between this script and the main one should leave the
26% pictures identical -- as long as the nodes stay uniform. Choose any other
27% node distribution and these weights are no longer the right ones, and the
28% "Poles" tab may well find real poles.
29% ---------------------------------------------------------------------------
31function w = bary_weights(x, d)
32n = numel(x) - 1;
33d = min(max(d, 0), n);
34w = zeros(1, n + 1);
35for k = 0:n
36 s = 0;
37 for i = max(0, k - d):min(k, n - d)
38 s = s + nchoosek(d, k - i);
39 end
40 w(k + 1) = (-1)^(k - d) * s;
41end
42end
44function r = bary_eval(x, y, w, t)
45sz = size(t);
46D = t(:) - x(:).';
47Q = w(:).' ./ D;
48r = (Q * y(:)) ./ sum(Q, 2);
49hit = find(any(D == 0, 2));
50for m = 1:numel(hit)
51 k = find(D(hit(m), :) == 0, 1);
52 r(hit(m)) = y(k);
53end
54r = reshape(r, sz);
55end