1function L = interpmat_1d(t,s)
2% INTERPMAT_1D interpolation matrix from nodes in 1D to any target nodes
3%
4% L = interpmat_1d(t,s) returns interpolation matrix taking values on nodes s
5% (a list of nodes) to target nodes t. It assumes smooth functions.
6% length(s) should be kept small, eg, 30 or less.
7%
8% Run without arguments does a self test (see test code for usage example).
9%
10% Notes: Computed in Helsing style, with centering and scaling for stability.
12% Barnett 7/17/16. Auto-centering & scaling for stability 12/23/21.
14if nargin==0, test_interpmat_1d; return; end
16cen = (max(s)+min(s))/2; hwid = (max(s)-min(s))/2; % affine s to [-1,1]
17s = (s-cen)/hwid;
18t = (t-cen)/hwid;
20p = numel(s); q = numel(t); s = s(:); t = t(:); % all col vecs
21n = p; % set the polynomial order we go up to
22V = ones(p,n); for j=2:n, V(:,j) = V(:,j-1).*s; end % polyval matrix on nodes
23R = ones(q,n); for j=2:n, R(:,j) = R(:,j-1).*t; end % polyval matrix on targs
24L = (V'\R')'; % backwards-stable way to do it (Helsing) See corners/interpdemo.m
26%%%%%%
27function test_interpmat_1d
28off = 4.3; % test centering and scaling
29sc = 1.7;
30x = sc*linspace(-1,1,16)' + off;
31f = @(x) sin(x + 0.7);
32data = f(x); % func on smooth (src) nodes
33t = sc*(2*rand(1000,1) - 1) + off; % cover same interval as the x lie
34uex = f(t); % col vec
35L = interpmat_1d(t,x);
36u = L * data;
37fprintf('max abs err for interp in [a,b] : %.3g\n',max(abs(u - uex)))
38fprintf('interp mat inf-norm = %.3g; max element size = %.3g\n',norm(L,inf),max(abs(L(:))))