Back to tutorial index

Splines and the geometry of curves


Topics in this lab

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

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

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)

Powered by MathJax

Published with MATLAB® 8.4