1% ---------------------------------------------------------------------------
2% Floater-Hormann barycentric rational interpolation
3% M. S. Floater and K. Hormann, Numer. Math. 107 (2007) 315-331
4%
5% This script is the method. Everything the four tabs draw comes out of the
6% functions below, so editing them changes the pictures.
7%
8% w = bary_weights(x, d) weights of equation (18) [required]
9% r = bary_eval(x, y, w, t) barycentric form, equation (1)[required]
10% [P, L] = local_blend(x, y, d, t) the blend of (4) and (5) [optional]
11% ---------------------------------------------------------------------------
13function w = bary_weights(x, d)
14% Equation (18). With J_k = { i : k-d <= i <= k } intersected with {0,...,n-d},
15%
16% w_k = sum_{i in J_k} (-1)^i prod_{j=i, j~=k}^{i+d} 1 / (x_k - x_j).
17%
18% Every k lies in the window of at most d+1 of the local polynomials, so this
19% costs O(n d^2) however the nodes are placed.
20n = numel(x) - 1;
21d = min(max(d, 0), n);
22w = zeros(1, n + 1);
23for k = 0:n
24 s = 0;
25 for i = max(0, k - d):min(k, n - d)
26 p = 1;
27 for j = i:(i + d)
28 if j ~= k
29 p = p / (x(k + 1) - x(j + 1));
30 end
31 end
32 s = s + (-1)^i * p;
33 end
34 w(k + 1) = s;
35end
36end
38function r = bary_eval(x, y, w, t)
39% Equation (1):
40%
41% r(t) = sum_k w_k f_k / (t - x_k) / sum_k w_k / (t - x_k).
42%
43% Berrut and Trefethen's advice is followed at the nodes themselves: if t is
44% exactly some x_k, return f_k rather than dividing by zero.
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
57function [P, L] = local_blend(x, y, d, t)
58% Equations (4) and (5), the construction the barycentric form above is a
59% rewriting of. P(i+1,:) is the polynomial p_i of degree at most d through the
60% d+1 points x_i, ..., x_{i+d}, and L(i+1,:) is the blending function lambda_i
61% normalised so that the columns of L sum to 1. Then r = sum_i L_i p_i.
62%
63% We build the blending functions from mu_i of equation (9),
64%
65% mu_i(t) = prod_{j<i} (t - x_j) * prod_{k>i+d} (x_k - t),
66%
67% rather than from lambda_i of equation (5) directly. The two differ by a
68% factor that does not depend on i, so they normalise to the same thing, but
69% mu_i is a polynomial and so has nothing to blow up at the nodes. Its sum is
70% the s(x) that Theorem 1 shows is positive, which is why L is well defined
71% everywhere.
72n = numel(x) - 1;
73d = min(max(d, 0), n);
74t = reshape(t, 1, []);
75m = n - d + 1;
77P = zeros(m, numel(t));
78Mu = zeros(m, numel(t));
79for i = 0:(n - d)
80 idx = (i + 1):(i + d + 1);
81 P(i + 1, :) = local_poly(x(idx), y(idx), t);
82 pr = ones(1, numel(t));
83 for j = 0:(i - 1)
84 pr = pr .* (t - x(j + 1));
85 end
86 for k = (i + d + 1):n
87 pr = pr .* (x(k + 1) - t);
88 end
89 Mu(i + 1, :) = pr;
90end
91L = Mu ./ sum(Mu, 1);
92end
94function p = local_poly(xi, yi, t)
95% Lagrange form of the polynomial of degree at most numel(xi)-1 through (xi, yi).
96p = zeros(1, numel(t));
97for k = 1:numel(xi)
98 b = ones(1, numel(t));
99 for j = 1:numel(xi)
100 if j ~= k
101 b = b .* (t - xi(j)) / (xi(k) - xi(j));
102 end
103 end
104 p = p + yi(k) * b;
105end
106end