5392320Interactive illustration of Floater-Hormann barycentric rational interpolationJeremy Magland 1function s = cubic_spline(x, y, dya, dyb, t)
2%CUBIC_SPLINE Clamped C^2 cubic spline interpolant, evaluated at t.
3%
4% This is the competitor in Tables 3 and 4 of the paper: a C^2 cubic spline
5% with clamped end conditions, i.e. with the first derivative at the two
6% end-points set to the corresponding derivative of f. Its error is O(h^4)
7% for f in C^4, the same order as the rational interpolant with d = 3.
8%
9% The moments M_i = s''(x_i) solve a tridiagonal system, which we solve with
10% the Thomas algorithm so that the cost stays O(n) even for the largest n in
11% the convergence study.
13x = x(:).';
14y = y(:).';
15n = numel(x) - 1;
16h = diff(x);
18lo = zeros(1, n + 1); % sub-diagonal
19di = zeros(1, n + 1); % diagonal
20up = zeros(1, n + 1); % super-diagonal
21rh = zeros(1, n + 1); % right-hand side
23% clamped end conditions
24di(1) = 2;
25up(1) = 1;
26rh(1) = 6 / h(1) * ((y(2) - y(1)) / h(1) - dya);
28for i = 2:n
29 hl = h(i - 1);
30 hr = h(i);
31 lo(i) = hl / (hl + hr);
32 di(i) = 2;
33 up(i) = hr / (hl + hr);
34 rh(i) = 6 * ((y(i + 1) - y(i)) / hr - (y(i) - y(i - 1)) / hl) / (hl + hr);
35end
37lo(n + 1) = 1;
38di(n + 1) = 2;
39rh(n + 1) = 6 / h(n) * (dyb - (y(n + 1) - y(n)) / h(n));
41% Thomas algorithm
42cp = zeros(1, n + 1);
43dp = zeros(1, n + 1);
44cp(1) = up(1) / di(1);
45dp(1) = rh(1) / di(1);
46for i = 2:n + 1
47 den = di(i) - lo(i) * cp(i - 1);
48 cp(i) = up(i) / den;
49 dp(i) = (rh(i) - lo(i) * dp(i - 1)) / den;
50end
51M = zeros(1, n + 1);
52M(n + 1) = dp(n + 1);
53for i = n:-1:1
54 M(i) = dp(i) - cp(i) * M(i + 1);
55end
57% evaluate: on [x_i, x_{i+1}] the spline is the usual cubic in the moments
58sz = size(t);
59t = t(:).';
60s = zeros(1, numel(t));
61for i = 1:n
62 if i == 1
63 m = t < x(2); % also catches t < x(1)
64 elseif i == n
65 m = t >= x(n); % also catches t > x(n+1)
66 else
67 m = (t >= x(i)) & (t < x(i + 1));
68 end
69 if ~any(m)
70 continue
71 end
72 tt = t(m);
73 hi = h(i);
74 ra = x(i + 1) - tt;
75 rb = tt - x(i);
76 s(m) = M(i) * ra.^3 / (6 * hi) + M(i + 1) * rb.^3 / (6 * hi) ...
77 + (y(i) - M(i) * hi^2 / 6) .* ra / hi ...
78 + (y(i + 1) - M(i + 1) * hi^2 / 6) .* rb / hi;
79end
80s = reshape(s, sz);
81end