Numerical Methods with Programming

Section 23.3 Using cubic splines

We do not go into the process of finding the coefficients of a cubic spline: it involves solving a linear system with only a few nonzero diagonals. In Matlab, the command spline computes a not-a-knot spline from given data, for example

x = 1:7;
y = [1 3 2 4 4 1 5];
s = spline(x, y);
disp(s.coefs)

Each row of the output contains the coefficients A, B, C, D of cubic polynomial (23.2.1) which is used on the interval \([x_{k}, x_{k+1}]\text{.}\) We could use these coefficients to evaluate the spline, but Matlab can do this automatically with ppval (meaning “Piecewise Polynomial VALues”). Specifically, ppval(s, t) evaluates the spline at a value (or vector of values) t. Thus,

t = linspace(min(x), max(x), 1000);
plot(t, ppval(s, t))

plots the spline computed above.

There is also an option to combine construction and evaluation of a spline into the same command spline, by feeding t as a third argument. That is, spline(x, y, t) constructs a spline from x, y and immediately evaluates it at t. For example:

t = linspace(min(x), max(x), 1000);
plot(t, spline(x, y, t))

We can also plot it with the original data points marked:

plot(t, spline(x, y, t), x, y, 'r*')

How else could we use the spline data? The formula (23.2.1) shows that \(S'(x_k) = C_k\text{,}\) which is the entry of the 3rd column of matrix s.coefs. Try s.coefs(:, 3)' and compare the values with the slope of the spline at the data points. This is another approach to numerical differentiation, different from what we did in Chapter 12: given discrete data, fit a spline to it and differentiate the spline.

Looking at the coefficients, we can also find the critical points of a spline, answering the question raised at the end of Section 23.1. Indeed, we have \(A, B, C, D\) coefficients of a cubic polynomial, so its derivative is a quadratic polynomial with coefficients \(3A, 2B, C\text{.}\) We can find its roots with roots but one must be careful to check if a root indeed falls into the subinterval on which this polynomial is used.

x = 1:7;
y = [1 3 2 4 4 1 5];

s = spline(x, y);
critpts = [];
for k=1:numel(x)-1
    r = roots([3 2 1].*s.coefs(k, 1:3));
    for j = 1:numel(r)
        if isreal(r(j)) && r(j) >= 0 && r(j) <= x(k+1)-x(k)
            critpts = [critpts, x(k)+r(j)];
        end
    end
end
t = linspace(min(x), max(x), 1000);
plot(t, ppval(s, t), x, y, 'r+', critpts, ppval(s, critpts), 'k*')