# Splines and the geometry of curves

## Topics in this lab

- Introduction
- Single polynomial interpolant (polyfit)
- Piecewise Polynomial (spline)
- First derivative (polyfit)
- First derivative (spline)
- Second derivative (polyfit)
- Second derivatives (spline)
- The Piecewise Cubic Hermite Interpolating Polynomial (PCHIP)
- First derivative (PCHIP)
- Second derivative (PCHIP)
- Computing tangents to the curve
- Lab 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+1 data points with a polynomial of degree N. 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+1); ydata = rand(1,N+1); 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); % Coefficients of polynomial of degree (N-1). ys = polyval(p,xs); hold on; plot(xs,ys,'b','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]);

Warning: Polynomial is badly conditioned. Add points with distinct X values, reduce the degree of the polynomial, or try centering and scaling as described in HELP POLYFIT.

Back to the top

## Piecewise Polynomial (spline)

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 spline is an
example of piecewise polynomial interpolant. We can simply use the
spline command with a single call. For example

ys = spline(xdata,ydata,xs); 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 the curve appears to have fewer oscillations than the single polynomial case.

As with the polyfit, we can also call the spline using a two step procedure. First, we compute the coefficients of the polynomial using the spline function, then we evaluate the polynomial using ppval.

pp_spline = spline(xdata,ydata)

pp_spline = form: 'pp' breaks: [1x14 double] coefs: [13x4 double] pieces: 13 order: 4 dim: 1

The pp_spline 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 piecewise polynomial, we use ppval

ys = ppval(pp_spline,xs);

This would produce the same y values that we computed above with the single call to spline.

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 (spline)

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

pp_spline.coefs

ans = 1.7189 -5.1881 3.8634 0.2290 1.7189 -1.2213 -1.0669 0.9133 -2.2680 2.7455 0.1056 0.1524 2.0906 -2.4884 0.3034 0.8258 -2.3459 2.3361 0.1862 0.5383 2.6330 -3.0775 -0.3841 0.9961 -2.3460 2.9987 -0.4447 0.0782 2.3944 -2.4152 0.0042 0.4427 -3.0751 3.1103 0.5389 0.1067 3.3069 -3.9862 -0.1349 0.9619 -2.3751 3.6452 -0.3972 0.0046 0.7988 -1.8358 0.9946 0.7749 0.7988 0.0076 -0.4117 0.8173

The rows of this array are the coefficients of the polynomial segments.
The form that each polynomial takes is given by
$$
p_k(x) = a_3(x-x_k)^3 + a_2(x-x_k)^2 + a_1(x-x_k) + a_0
$$
where the $x_k$ are the breakpoints
typically defined by the x data points.
We then evaluate the interpolant using ppval.

Using the structure pp_spline in combination with
polyder, we can compute the coefficients of the
of each the derivative of each
polynomial segment. Using these coefficients, we can then create a new
piecewise polynomial structure which we can easily evaluate.

% loop over each interval to get the derivatives in that interval. coeffs_deriv = zeros(N,3); % store the coefficients of the derivative for k = 1:N, c = pp_spline.coefs(k,:); % pp is a structure coeffs_deriv(k,:) = polyder(c); end pp_deriv = mkpp(xdata,coeffs_deriv); ys_deriv = ppval(pp_deriv,xs); ydata_deriv = ppval(pp_deriv,xdata); 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 derivatives (spline)

The cubic spline has a continuous second derivative.

% loop over each interval to get the derivatives in that interval. coeffs_deriv2 = zeros(N,2); for k = 1:N, c = pp_deriv.coefs(k,:); % pp is a structure coeffs_deriv2(k,:) = polyder(c); end pp_deriv2 = mkpp(xdata,coeffs_deriv2); ys_deriv2 = ppval(pp_deriv2,xs); ydata_deriv2 = ppval(pp_deriv2,xdata); 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

## The Piecewise Cubic Hermite Interpolating Polynomial (PCHIP)

While both the single polynomial and the spline function produce smooth curves that interpolate our data, we may not be so happy with the fact that in some cases, the curve may introduce new extrema that were not present in the original data. This is especially problematic if the data is only physically meaningful in a limited range, say, [0,1]. In this situation, Matlab provides a useful command that prevents the introduction of new extrema into the data. This is called the pchip and can be used in a manner very similar to the spline comand. Below, we illustrate the pchip interpolant and the first two derivatives. You'll notice that the second derivative is no longer continuous.

ys = pchip(xdata,ydata,xs); clf; plot(xs,ys,'g','linewidth',2); hold on; plot(xdata,ydata,'k.','markersize',30); plot(xlim,[0 0],'k--'); plot(xlim,[1 1],'k--'); title('Interplating function (pchip)','fontsize',18); ylim([-1 2]);

Notice that all of the values are between 0 and 1. The values returned from evaluating a pchip curve will never exceed range of the original data.

Back to the top

## First derivative (PCHIP)

The pchip spline, unlike the spline derivative,
uses derivatives that are computed strictly from the data. Derivative
values at extrema are set to 0, whereas derivatives at intermediate data
points are taken to be a weighted average of the two one-sided slopes
based on y-data.

First, we have to create a spline for the first and second derivatives.

pp_pchip = pchip(xdata,ydata); coeffs_deriv = zeros(N,3); coeffs_deriv2 = zeros(N,2); for k = 1:N, c = pp_pchip.coefs(k,:); coeffs_deriv(k,:) = polyder(c); coeffs_deriv2(k,:) = polyder(polyder(c)); end pp_deriv = mkpp(xdata,coeffs_deriv); pp_deriv2 = mkpp(xdata,coeffs_deriv2);

Then we can evaluate and plot these splines representing the derivatives.

clf; ys_deriv = ppval(pp_deriv,xs); ydata_deriv = ppval(pp_deriv,xdata); plot(xs,ys_deriv,'g','linewidth',2); hold on; plot(xdata,ydata_deriv,'k.','markersize',30); plot(xlim,[0 0],'k--'); title('First derivative (pchip)','fontsize',18);

Back to the top

## Second derivative (PCHIP)

The second derivative of the PCHIP spline is not continuous, as you can see from this graph

clf; ys_deriv2 = ppval(pp_deriv2,xs); ydata_deriv2 = ppval(pp_deriv2,xdata); plot(xs,ys_deriv2,'g.','markersize',10); hold on; plot(xdata,ydata_deriv2,'k.','markersize',30); title('Second derivative (pchip)','fontsize',18);

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.

f_spline = spline(xdata,ydata); coeffs_deriv = zeros(N,3); coeffs_deriv2 = zeros(N,2); for k = 1:N, c = f_spline.coefs(k,:); coeffs_deriv(k,:) = polyder(c); coeffs_deriv2(k,:) = polyder(polyder(c)); end f_deriv = mkpp(xdata,coeffs_deriv); f_deriv2 = mkpp(xdata,coeffs_deriv2); ys = ppval(f_spline,xs); clf; plot(xs,ys,'r','linewidth',2); hold on; plot(xdata,ydata,'k.','markersize',30);

Now we will use the formulas to compute the tangents to the curve at our data points.

```
% Compute tangent vector
r = [xdata; ydata];
rp = [ones(1,N+1); ppval(f_deriv,xdata)];
g = sqrt(dot(rp,rp));
T = bsxfun(@times,rp,1./g);
```

We can plot the tangents as arrows using the quiver command.

quiver(xdata,ydata,T(1,:),T(2,:),0.25,'k','linewidth',2); hold on; plot(xdata,ydata,'k.','markersize',30); title('Tangents (spline)','fontsize',18); xlim([-5 5]); ylim([-.5 1.5]);

Back to the top

## Lab 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.4