1function [r, poles] = classical_rational(x, y, t)
2%CLASSICAL_RATIONAL "Classical" rational interpolation p_M / q_N, M + N = n.
3%
4% This is the construction the paper's introduction describes and rejects:
5% fit the values f(x_i) with a quotient of polynomials of degrees M and N
6% with M + N = n, taking M = N = n/2 when n is even. It is the method with
7% "no control over the occurrence of poles in the interval of interpolation",
8% and this routine returns those poles so that they can be drawn.
9%
10% The interpolation conditions p(x_i) - y_i q(x_i) = 0 are linear in the
11% coefficients, so the coefficient vector is a null vector of an
12% (n+1) x (n+2) matrix, which we take from the last right singular vector.
13% Everything is done in a variable rescaled to [-1, 1], since the monomial
14% basis on the original interval is badly conditioned.
15%
16% Note that solving the linearised conditions does not guarantee that the
17% quotient actually interpolates: a common root of p and q at some x_i (an
18% "unattainable point") is possible. That is a further wrinkle of the
19% classical method, not a bug here.
21x = x(:).';
22y = y(:).';
23n = numel(x) - 1;
24M = ceil(n / 2);
25N = n - M;
27a = min(x);
28b = max(x);
29c0 = (a + b) / 2;
30sc = (b - a) / 2;
31xs = (x - c0) / sc;
33V = zeros(n + 1, M + N + 2);
34for i = 1:(n + 1)
35 V(i, 1:(M + 1)) = xs(i).^(M:-1:0);
36 V(i, (M + 2):end) = -y(i) * xs(i).^(N:-1:0);
37end
39[~, ~, W] = svd(V);
40c = W(:, end).';
41pc = c(1:(M + 1));
42qc = c((M + 2):end);
44ts = (t - c0) / sc;
45r = polyval(pc, ts) ./ polyval(qc, ts);
46r = reshape(r, size(t));
48% real poles = real roots of q, mapped back to the original variable
49poles = [];
50tol = 1e-13 * max(abs(qc));
51k0 = find(abs(qc) > tol, 1);
52if ~isempty(k0) && numel(qc) - k0 >= 1
53 z = roots(qc(k0:end));
54 z = z(:).';
55 keep = abs(imag(z)) < 1e-7 * max(1, max(abs(z)));
56 poles = sort(real(z(keep)) * sc + c0);
57end
58end