5392320Interactive illustration of Floater-Hormann barycentric rational interpolationJeremy Magland 1# barycentric-rational
3An interactive illustration of
5> M. S. Floater and K. Hormann, **Barycentric rational interpolation with no poles
6> and high rates of approximation**, *Numerische Mathematik* **107** (2007) 315–331.
7> [doi:10.1007/s00211-007-0093-y](https://doi.org/10.1007/s00211-007-0093-y)
9The method itself is a MATLAB-syntax script you can edit in the page. It runs in
10your browser through [numbl](https://numbl.org); there is no server and nothing
11to install.
13## What the paper says
15Interpolating a function at a given set of points is easy to do badly. The
16degree-*n* polynomial through *n*+1 equally spaced points diverges as *n* grows,
17which is Runge's example. Classical rational interpolation, fitting a quotient
18p<sub>M</sub>/q<sub>N</sub> with M + N = n, often approximates better but offers
19no control over where the poles land, and they land inside the interval.
21Floater and Hormann's construction is short. Fix an integer *d* with 0 ≤ *d* ≤ *n*.
22For each *i* let p<sub>i</sub> be the polynomial of degree at most *d* through the
23*d*+1 points x<sub>i</sub>, …, x<sub>i+d</sub>, and blend those local polynomials
24together:
26$$r(x) = \frac{\sum_{i=0}^{n-d} \lambda_i(x)\, p_i(x)}{\sum_{i=0}^{n-d} \lambda_i(x)},
27\qquad \lambda_i(x) = \frac{(-1)^i}{(x - x_i)\cdots(x - x_{i+d})}.$$
29The results are that *r* has **no poles anywhere on the real line** for any *d*
30(Theorem 1), that its error is **O(h<sup>d+1</sup>)** for *d* ≥ 1 **whatever the
31node distribution**, as long as *f* is smooth enough (Theorem 2), and that *r* can
32be rewritten in the barycentric form
34$$r(x) = \sum_{k=0}^{n} \frac{w_k}{x - x_k} f(x_k) \Big/ \sum_{k=0}^{n} \frac{w_k}{x - x_k}$$
36with explicit weights (equation 18), which is cheap to evaluate. On a uniform mesh
37those weights are integers, nearly all equal, differing only near the two ends:
381, 4, 7, 8, …, 8, 7, 4, 1 for *d* = 3. That small change at the ends is what lifts
39the approximation order from O(h) to O(h<sup>4</sup>).
41The *d* = 0 case is Berrut's earlier interpolant, and *d* = *n* is ordinary
42polynomial interpolation, so the family interpolates between the two.
44## What the page shows
46Four tabs, all driven by the same script:
48- **Interpolant** — *f*, the rational interpolant *r*, and the nodes, with the
49 degree-*n* polynomial and a clamped C² cubic spline as optional overlays, plus
50 the pointwise error underneath.
51- **Blending & weights** — the *n*−*d*+1 local polynomials and the normalised
52 blending functions λ<sub>i</sub>, which sum to 1 everywhere but have oscillating
53 tails and no local support; and a stem plot of the barycentric weights, with the
54 integers of Section 4 read off when the mesh is uniform.
55- **Poles** — the denominator of *r* on the real line, drawn as a signed *n*-th
56 root so that its sign and zeros survive the enormous dynamic range; all of its
57 roots plotted in the complex plane, none of them touching the real axis; and the
58 classical p<sub>M</sub>/q<sub>N</sub> alongside, with the poles it does put in
59 the interval.
60- **Convergence** — max error against *n* on log-log axes for a range of *d*, with
61 the measured orders tabulated. With Runge's function on uniform nodes this
62 reproduces Table 1 of the paper, and with the spline turned on, Tables 3 and 4.
64## The script
66The editor holds the whole method. The app only asks it for two functions, and
67optionally a third:
69```matlab
70w = bary_weights(x, d) % the weights, equation (18)
71r = bary_eval(x, y, w, t) % the barycentric form, equation (1)
72[P, L] = local_blend(x, y, d, t) % the blend of (4) and (5) [optional]
73```
75Anything that satisfies that contract will drive all four tabs, which is the point
76of the alternative scripts in the **method** dropdown:
78| script | what it does |
79|---|---|
80| Floater-Hormann | the paper |
81| Berrut (d = 0) | weights (−1)<sup>k</sup>; set the nodes to **paired** to see why Theorem 3 needs a bounded mesh ratio |
82| Integer weights | Section 4's closed form; identical to the first script while the mesh stays uniform, and not otherwise |
83| Polynomial (d = n) | the Lagrange weights of equation (2); no poles, but Runge divergence |
84| Equal weights | drop the alternating signs and a pole appears in every interval |
85| Random weights | the generic barycentric rational interpolant: interpolates, has poles |
87The last two are the counterpart to Schneider and Werner's theorem, quoted in the
88paper, that a pole-free barycentric rational interpolant must have weights that
89alternate in sign.
91## Running it locally
93```bash
94npm install
95npm run dev # http://localhost:5173
96npm run build # type-check and bundle to dist/
97```
99Two test suites, neither of which needs a browser to be watched:
101```bash
102npm run test:matlab # the .m layer against the paper's tables, via the numbl CLI
103npm run test:matlab -- --full # also n = 640
104npm run test:browser # the built app in headless Chrome
105```
107`test:matlab` runs the MATLAB layer outside the browser through a clone of
108[numbl](https://github.com/flatironinstitute/numbl) (set `NUMBL_DIR`; it defaults
109to `~/src/numbl`) and checks it against the published numbers: the integer weight
110patterns of Section 4 for *d* = 0…4, the error columns of Tables 1, 3 and 4, the
111absence of real roots for every *d* across five node distributions, and the
112identity r = Σ L<sub>i</sub> p<sub>i</sub> relating equations (1) and (4).
114One entry of the paper's Table 1 does not reproduce: the sine row at *n* = 20 is
115printed as 3.9e−05, but the order 5.5 printed beside it implies 1.7e−2 / 2<sup>5.5</sup>
116= 3.8e−04, which is what we get, and every other entry in the row matches to two
117figures. We take it as a misprint.
119## How it is put together
121- `src/matlab/driver.m` — reads `params.json`, calls the method script, writes
122 `out.json`. It is prepended to whatever is in the editor, so the script's
123 functions become its local functions.
124- `src/matlab/lib/*.m` — the parts that are not the method: node distributions,
125 test functions, the cubic spline, the classical rational interpolant.
126- `src/methods/*.m` — the scripts in the dropdown.
127- `src/engine/` — the numbl session. One session per script; parameter changes
128 only rewrite `params.json` and re-run, which is fast enough to drive a slider.
129- `src/plot/` — a small SVG plotting layer. The palette is documented and was
130 checked with a colour-vision validator rather than by eye.
131- `src/panels/` — one component per tab.
133## License
135Apache-2.0, matching numbl.