# Splines and the geometry of curves

## Introduction

Cubic splines are piecewise polynomials which are continuous and have continuous derivatives. In this lab, we will investigate how to use the Matlab functions pchip and spline. As we will see, pchip is good for interpolating physical data which should be constrained in some way (should always be positive, or between 0 and 1, for example) while the spline function is ideally suited for geometric data.

clear all;
close all;

format short;

Back to the top

## Single polynomial interpolant (polyfit)

Before starting with splines, we illustrate how to fit a single polynomial through a set of data points. Recall that we can uniquely interpolate N data points with a polynomial of degree N-1. As an example, suppose we wish to put a smooth curve through a set of random data points.

N = 13;
xdata = linspace(-5,5,N);
ydata = rand(1,N);

xs = linspace(-5,5,500);

plot(xdata,ydata,'k.','markersize',30);
hold on;
plot(xdata,0*xdata,'k--');
plot(xdata,0*xdata + 1,'k--');

ylim([-1 2]);

We first construct the polynomial using the polyfit function. This function returns the coefficients of the interpolating polynomial. We then evaluate the polynomial using polyval.

p = polyfit(xdata,ydata,N-1);   % Coefficients of polynomial of degree (N-1).
ys = polyval(p,xs);

hold on;
plot(xs,ys,'r','linewidth',2);
plot(xdata,ydata,'.k','markersize',30);
plot(xdata,0*xdata,'k--');
plot(xdata,0*xdata + 1,'k--');
title('Interpolating function (polyfit)','fontsize',18);
ylim([-1 2]);

Notice that even though the original data is between 0 and 1, the interpolant has values that are negative, or wildly exceed 1.

Back to the top

## Piecewise Cubic Hermite Interpolating Polynomial (pchip)

We can improve the behavior of the single polynomial interpolant by adjusting the location of the points in xdata. But often, we do not have a choice as to where to place the data points. A second option is to construct a piecewise polynomial. A common approach is to put a cubic polynomial over each interval defined by the N data points. The function pchip is one example of piecewise polynomial interpolant, the Piecewise Cubic Hermite Interpolating Polynomial. As with the polyfit, we use a two step procedure. First, we compute the coefficients of the polynomial using the pchip function, then we evaluate the polynomial using ppval.

p_pchip = pchip(xdata,ydata)
p_pchip =

form: 'pp'
breaks: [1x13 double]
coefs: [12x4 double]
pieces: 12
order: 4
dim: 1



The p_pchip is a Matlab structure that stores several pieces of information. In this case, the structure stores the break points (defined by xdata), the coefficients of each polynomial in the segments defined by the breakpoints, the dimension of the polynomial, the number of pieces, and the order.

To evaluate the polynomial, we use ppval

ys = ppval(p_pchip,xs);

Finally, we plot the resulting polynomial.

% Plot the resulting polynomial
clf;
plot(xs,ys,'r','linewidth',2);
hold on;
plot(xdata,ydata,'k.','markersize',30);
plot(xdata,0*xdata,'k--');
plot(xdata,0*xdata + 1,'k--');
title('Interpolating function (pchip)','fontsize',18);
ylim([-1 2]);

Notice that using the PCHIP interpolant, the values of the interpolant are all between 0 and 1.

Back to the top

## Cubic spline (spline)

Another type of cubic piecewise polynomial can be constructed using the spline function. The approach to using the spline is the same as for using pchip.

p_spline = spline(xdata,ydata);
ys = ppval(p_spline,xs);

% Plot the polynomial
clf;
plot(xs,ys,'r','linewidth',2);
hold on;
plot(xdata,ydata,'k.','markersize',30);
plot(xdata,0*xdata,'k--');
plot(xdata,0*xdata + 1,'k--');
title('Interpolating function (spline)','fontsize',18);
ylim([-1 2]);

Notice that spline is smoother, but that the it may exceed the range (0,1).

Back to the top

## First derivative (polyfit)

To compute derivatives of the polynomial interpolant, we need to be able to compute derivatives of a polynomial described by a set of coefficients. For example, to compute the derivative of the polynomial obtained using polyfit, we use polyder as

clf;
p = polyfit(xdata,ydata,N-1);

p_deriv = polyder(p);
ys_deriv = polyval(p_deriv,xs);
ydata_deriv = polyval(p_deriv,xdata);

plot(xs,ys_deriv,'b','linewidth',2);
hold on;
plot(xdata,ydata_deriv,'k.','markersize',30);
title('First derivative (polyfit)','fontsize',18);
ylim([-10 10]);

Notice the scale on the y-axis. The derivative is quite large.

Back to the top

## First derivative (pchip)

To compute the derivative of the piecewise polynomial, we first need access to the coefficients of each piecewise segment. These coefficients can be accessed from the structure p_pchip.

p_pchip.coefs
ans =

0.1574   -0.3934    0.3278    0.8147
2.6916   -3.3644         0    0.9058
-2.7178    3.3972         0    0.1270
0.3345   -0.6834         0    0.9134
1.2117   -1.2493   -0.4421    0.6324
-0.2519    0.4705         0    0.0975
0.0069    0.0694    0.2594    0.2785
-0.8331    0.8182    0.3895    0.5469
-0.0005   -0.0099    0.0174    0.9575
2.7899   -3.4874         0    0.9649
-2.8097    3.5121         0    0.1576
-0.0232    0.0000         0    0.9706



The rows of this array are the coefficients of the polynomial segments. The form that each polynomial takes is given by $$p_j(x) = c_3(x-b_j)^3 + c_2(x-b_j)^2 + c_1(x-b_j) + c_0$$ where the $b_j$ are the breakpoints typically defined by the x data points. We then evaluate the interpolant using ppval.

Using the structure p_pchip in combination with polyder, we can compute the derivative of each polynomial segment. Below, we plot the derivative of the PCHIP polynomial.

ys_deriv = 0*xs;

% loop over each interval to get the derivatives in that interval.
for j = 1:N-1,
c = p_pchip.coefs(j,:);   % pp is a structure

m = xdata(j) < xs & xs <= xdata(j+1);
c_deriv = polyder(c);
ys_deriv(m) = polyval(c_deriv,xs(m)-xdata(j));
ydata_deriv(j) = polyval(c_deriv,0);
end

clf;
plot(xs,ys_deriv,'r','linewidth',2);
hold on;
plot(xdata,ydata_deriv,'k.','markersize',30);
plot(xdata,0*xdata,'k--');
title('First derivative (pchip)','fontsize',18);
ylim([-2 2]);

Notice that the derivatives at several of the data points (or knots) is exactly zero. This is because pchip approximates derivatives at the data points based on an alogrithm which assigns a derivative of 0 to data points which represent a local maximum or minimum. See Numerical Computing with Matlab for a complete discussion of the pchip algorithm.

Back to the top

## First derivative (spline)

We can compute the piecewise derivatives of the spline function in the same way.

p_spline = spline(xdata,ydata);
ys_deriv = 0*xs;

% loop over each interval to get the derivatives in that interval.
for j = 1:N-1,
c = p_spline.coefs(j,:);   % pp is a structure

m = xdata(j) < xs & xs <= xdata(j+1);
c_deriv = polyder(c);
ys_deriv(m) = polyval(c_deriv,xs(m)-xdata(j));
ydata_deriv(j) = polyval(c_deriv,0);
end

clf;
plot(xs,ys_deriv,'r','linewidth',2);
hold on;
plot(xdata,ydata_deriv,'k.','markersize',30);
plot(xdata,0*xdata,'k--');
title('First derivative (spline)','fontsize',18);
ylim([-2 2]);

Back to the top

## Second derivative (polyfit)

Continuing in this fashion, we can plot the second derivatives of each of these interpolants.

clf;
p = polyfit(xdata,ydata,N-1);

p_deriv2 = polyder(polyder(p));
ys_deriv2 = polyval(p_deriv2,xs);
ydata_deriv2 = polyval(p_deriv2,xdata);

plot(xs,ys_deriv2,'b','linewidth',2);
hold on;
plot(xdata,ydata_deriv2,'k.','markersize',30);
title('Second derivative (polyfit)','fontsize',18);
ylim([-10 10]);

Back to the top

## Second derivative (pchip)

With pchip, it is obvious that the second derivative of the piecewise interpolating polynomial is not continuous.

pp = pchip(xdata,ydata);
ys_deriv2 = 0*xs;

% loop over each interval to get the derivatives in that interval.
for j = 1:N-1,
c = pp.coefs(j,:);   % pp is a structure

m = xdata(j) < xs & xs <= xdata(j+1);
c_deriv2 = polyder(polyder(c));
ys_deriv2(m) = polyval(c_deriv2,xs(m)-xdata(j));
ydata_deriv2(j) = polyval(c_deriv2,0);
end

clf;
plot(xs,ys_deriv2,'r.','markersize',10);
hold on;
plot(xdata,ydata_deriv2,'k.','markersize',30);
title('Second derivative (pchip)','fontsize',18);
ylim([-10 10]);

Back to the top

## Second derivatives (spline)

The cubic spline, however, has a continuous second derivative.

pp = spline(xdata,ydata);
ys_deriv2 = 0*xs;

% loop over each interval to get the derivatives in that interval.
for j = 1:N-1,
c = pp.coefs(j,:);   % pp is a structure

m = xdata(j) < xs & xs <= xdata(j+1);
c_deriv2 = polyder(polyder(c));
ys_deriv2(m) = polyval(c_deriv2,xs(m)-xdata(j));
ydata_deriv2(j) = polyval(c_deriv2,0);
end

clf;
plot(xs,ys_deriv2,'r.','markersize',10);
hold on;
plot(xdata,ydata_deriv2,'k.','markersize',30);
title('Second derivative (spline)','fontsize',18);
ylim([-10 10]);

Back to the top

## Computing tangents to the curve

We can compute tangents to the curve using a formula for the tangent. For a function, we assume the data points have the form ${\bf r}(x) = (x,f(x))$. Then we can apply the formula for the tangent, given by $${\bf T}(x) = \frac{{\bf r}'(x)}{g(x)}$$ where $$g(x) = \parallel \!\! {\bf r}'(x) \!\! \parallel = \sqrt{1 + f'(x)^2}$$ We will test this out on the interpolating spline.

clf;
pp = spline(xdata,ydata);
ys = ppval(pp,xs);

plot(xs,ys,'r.','markersize',10);
hold on;

% Compute the Tangents at the knots (xdata,ydata)
T = zeros(2,N+1);
for j = 1:N-1,
c = pp.coefs(j,:);   % pp is a structure

c_deriv = polyder(c);
fp = polyval(c_deriv,0);
rp = [1; fp];
g = norm(rp,2);  % = sqrt(1 + f'(x)^2)
T = rp/g;
plot([xdata(j) xdata(j)+T(1)],[ydata(j) ydata(j)+T(2)],'k','linewidth',2);
end

plot(xdata,ydata,'k.','markersize',30);
title('Tangents (spline)','fontsize',18);
xlim([-5 5]);
ylim([-1.5 1.5]);
% daspect([1 1 1]);

Back to the top

## Exercises

1. Use the functions polyval and polyder to plot the following polynomial and its derivative. $$p(x) = 5x^3 - 2x^2 + 4x - 10$$
2. Use the functions polyval and polyder to plot the following polynomial and its derivative. $$p(x) = 5(x-1)^3 - 2(x-1)^2 + 4(x-1) - 10$$
3. Using the norm function to compute the 2-norm $\parallel \cdot \parallel$, show that $$\parallel \!\! {\bf v} \! \! \parallel = \sqrt{v_1^2 + v_2^2 + v_3^2}$$ where ${\bf v} = (-1,2,5)$. You can compute the 2-norm of a vector (i.e. the length of the vector) with the command norm(v,2).
4. Construct a cubic spline through the polynomial $$p(x) = x(x^2-1)$$ at the points $[-2,-1,0,1,2]$. Use the function ppval to interpolate the spline at a set of points on the interval $[-2,2]$. Use polyder to compute the derivative of the polynomial. You can start with this code fragment.
xdata = [-2,-1,0,1,2];
ydata = xdata.*(xdata.^2 - 1);
pp = spline(xdata,ydata);
xs = linspace(-2,2,100);
ys = ppval(pp,xs);


Be sure you understand the structure pp. Plot the original polynomial and its derivative on top of your spline results to see that your spline results are correct.
5. Construct and plot cubic spline through data points given by
N = 8;
xdata = linspace(-5,5,N);
ydata = -2 + 4*rand(1,N);

6. Using the the data points from the previous example, compute and plot the tangents and normals to the curves. Use the formulas $${\bf T}(x) = \frac{{\bf r}'(x)}{g(x)}$$ and $${\bf N}(x) = \frac{g(x) {\bf r}''(x) - g'(x){\bf r}'(x)}{g(x)^3 \kappa(x)}$$ where $$g(x) = \parallel \! \! {\bf r}'(x)\! \! \parallel \quad \mbox{and} \quad \kappa(x) = \frac{\parallel \!\! {\bf r}''(x) \times {\bf r}'(x) \!\! \parallel}{g(x)^3}$$ To compute the 2-norm $\parallel {\bf v} \parallel$, use norm(v,2). To compute the cross product ${\bf u}\times {\bf v}$, use the command cross(u,v).
Plot the tangents and the normals. Your final plot should look like :

7. Use the dot function to show that the normals and the tangent vectors you created are perpendicular. That is, show that $${\bf T} \cdot {\bf N} = 0$$.
8. Show that the binormal vector, computed as $${\bf B} = {\bf T} \times {\bf N}$$ has zeros in the x- and y- components, is only non-zero in the z-component.
9. Show that $${\bf T} \cdot {\bf B} = {\bf N} \cdot {\bf B} = 0$$
10. Show that $$\parallel \!\! {\bf B} \! \! \parallel = \parallel \!\! {\bf N} \! \! \parallel = \parallel \!\! {\bf T} \! \! \parallel = 1$$