1function x = nodes_of(kind, a, b, n, seed)
2%NODES_OF Build n+1 interpolation nodes a = x_0 < x_1 < ... < x_n = b.
3%
4% Floater and Hormann's Theorem 2 gives the rate O(h^(d+1)) for d >= 1
5% regardless of how the nodes are distributed, so it is worth being able to
6% distribute them badly. The 'paired' family below does exactly that: it
7% pulls every second node close to its left neighbour, which drives the local
8% mesh ratio beta of Theorem 3 up and makes the d = 0 case (Berrut's
9% interpolant) misbehave while d >= 1 carries on unaffected.
11k = 0:n;
12switch kind
13 case 'uniform'
14 x = a + (b - a) * k / n;
15 case 'chebyshev'
16 % Chebyshev-Gauss-Lobatto points, clustered at both ends
17 x = (a + b) / 2 - (b - a) / 2 * cos(pi * k / n);
18 case 'random'
19 rng(seed);
20 u = sort(rand(1, max(0, n - 1)));
21 x = [a, a + (b - a) * u, b];
22 case 'paired'
23 x = a + (b - a) * k / n;
24 h = (b - a) / n;
25 % indices 2,4,... are the nodes x_1, x_3, ... (1-based vs 0-based)
26 x(2:2:end) = x(2:2:end) - 0.9 * h;
27 case 'graded'
28 % quadratically graded, clustered at the left end
29 x = a + (b - a) * (k / n).^2;
30 otherwise
31 error('nodes_of: unknown node distribution ''%s''', kind);
32end
33x = x(:).';
34end