/ concept-collection / barycentric-rational
Sign in
concept-collection / barycentric-rational
barycentric-rational / src / methods / random.m
35 lines · 1.3 KBCodeBlameHistory
5392320Interactive illustration of Floater-Hormann barycentric rational interpolationJeremy Magland 1% ---------------------------------------------------------------------------
2% Random weights: the generic barycentric rational interpolant
3%
4% Berrut and Mittelmann's observation, quoted in the paper's introduction, is
5% that *every* rational interpolant whose numerator and denominator have degree
6% at most n can be written in the barycentric form (1) for some real weights
7% w_0, ..., w_n. So the whole difficulty of the subject is choosing them. This
8% script chooses them at random.
9%
10% The result still interpolates -- that is free -- but the denominator now
11% changes sign wherever it pleases, and the Poles tab finds real poles inside
12% the interval. Reroll with the d slider, which is used here only to reseed.
14% This is the situation the paper's construction is a way out of: weights that
15% are known in advance to give a pole-free interpolant, and one whose
16% approximation order can be raised at will.
17% ---------------------------------------------------------------------------
19function w = bary_weights(x, d)
20rng(1 + d);
21w = randn(1, numel(x));
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
moveopenescclose