# Splines and the geometry of curves

## Topics in this lab

- Introduction
- Single polynomial interpolant (polyfit)
- Piecewise Cubic Hermite Interpolating Polynomial (pchip)
- Cubic spline (spline)
- First derivative (polyfit)
- First derivative (pchip)
- First derivative (spline)
- Second derivative (polyfit)
- Second derivative (pchip)
- Second derivatives (spline)
- Computing tangents to the curve
- Exercises

## 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

- Use the functions polyval and polyder to plot the following polynomial and its derivative. $$p(x) = 5x^3 - 2x^2 + 4x - 10$$
- 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 $$
- 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).
- 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. - Construct and plot cubic spline through data points given by
N = 8; xdata = linspace(-5,5,N); ydata = -2 + 4*rand(1,N);

- 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 :

- 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$$.
- 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. - Show that $${\bf T} \cdot {\bf B} = {\bf N} \cdot {\bf B} = 0$$
- Show that $$ \parallel \!\! {\bf B} \! \! \parallel = \parallel \!\! {\bf N} \! \! \parallel = \parallel \!\! {\bf T} \! \! \parallel = 1 $$

Compare your answers with the solutions.

Back to the top

## Get the code

Do you want to try the above code fragments on your own? Download the Matlab script that produces this page here. (lab_15.m)

Published with MATLAB® 8.2