Interpolation in MATLAB

I'd like to introduce a new guest blogger - John D'Errico - an applied mathematician, now retired from Eastman Kodak, where he used MATLAB for over 20 years. Since then, MATLAB is still in his blood, so you will often find him answering questions on the newsgroup and writing new utilities to add to MATLAB Central.

Contents

I'll assume you have some data points through which you wish to pass a curve, interpolating your data. (Initially, I will only talk about problems with one independent variable.) In these coming blogs, I'll try to show some ways to do exactly this, i.e., find a curve that passes through your data. Along the way I'll try to give some pointers on curve fitting, interpolation, modeling, approximation, etc.

Polynomial regression

A valid question for some to ask is why start out with a discussion about polynomial regression , when we really wanted to talk about interpolation. Many people mistake the ideas of interpolation with the approximation produced by a regression model, calling both of these things interpolation. So I'm starting out with some discussion about what interpolation is not. Plus, I want to assure an understanding of polynomials, since many of the tools for interpolation are polynomial based in some way.

Let us start by creating some data. An exponential is a good place to start, a simple curve shape that is easy to fit.

x = -1:.1:1;
y = exp(x);

Plot your data

It is always a good idea to plot your data. In fact, I'll suggest that you should plot everything. Plots are useful, since your eye and brain are splendid at things like pattern recognition. Only you know your data, as the scientist, engineer, or analyst. You will always benefit if you can employ your knowledge of a system as part of the modeling process.

Next, always think in advance about your goals for any model.

For example, I once had a problem where I knew I had some significant noise in our process, but I chose not to do any smoothing anyway. Any such smoothing would also have smoothed out some potentially important features of our process. Since I could survive with the noise in my interpolant, I chose the lesser evil.

Always know your goals for any such task.

plot(x,y,'bo')
xlabel X
ylabel Y
grid on
title 'Exponential data'

This is a nice, well-behaved function. It is of the form y = f(x). I'll pretend for the moment that I have no idea what was in the underlying functional relationship.

One thing I learned in some long past calculus course was that a Taylor series will provide an approximation for many functions. Polynomials are useful things. They are simple to use, simple to build, simple to work with. And a truncated Taylor series is basically a nice, simple polynomial.

So we will start with a linear polynomial approximation for this curve, built using polyfit . This is a utility provided in MATLAB to estimate a polynomial model using linear regression techniques. We could also use many other tools to build our polynomial model, but polyfit is a useful one, and easy to use.

A linear, or first degree polynomial (many use the words "order" and "degree" interchangeably), might be written mathematically as y(x) = a1*x + a2. In MATLAB we will merely store the coefficients, as a vector [a1,a0]. Note that a polynomial in MATLAB has it's coefficients stored with the highest order term first.

P1 = polyfit(x,y,1)

P1 =
    1.1140    1.1937

We can evaluate the polynomial with polyval.

yhat = polyval(P1,x);

plot(x,y,'bo',x,yhat,'r-')
xlabel X
ylabel Y
grid on
title 'Linear polynomial fit'

Look at the residuals

In fact, I'll claim the relationship we are modeling is not terribly well represented by a linear model. Depending on your needs for this model, you might have decided differently.

When you build a regression model, look at the residuals.

ALWAYS PLOT EVERYTHING!

Plot your residuals. In general, some good ways to plot the residuals are versus

Look at your plots. Think about what you see there. Compare it to your expectations.

Since this data was very simply generated, I'll dispense with some of those plots for brevity. Note that the residuals for this linear fit look vaguely like a quadratic polynomial. This is often the case when there is lack of fit in a polynomial. That lack of fit often looks like the first term we truncated from the Taylor series.

res = y - yhat;
plot(x,res,'bo')
xlabel X
ylabel Residuals
grid on
title 'Residuals for the linear fit'

If the residuals looked vaguely parabolic in shape, then it might make sense to use a second order (quadratic) polynomial for our fit.

P2 = polyfit(x,y,2)
yhat = polyval(P2,x);
plot(x,y,'bo',x,yhat,'r-')
xlabel X
ylabel Y
grid on
title 'Quadratic polynomial fit'

P2 =
    0.5402    1.1140    0.9956

Note that the residuals here look vaguely like a cubic polynomial, although they are much smaller in magnitude than the previous fit. Again, at each step as we increase the order of the model, the residuals will often tend to look much like a polynomial of the next higher order.

res = y - yhat;
plot(x,res,'bo')
xlabel X
ylabel Residuals
grid on
title 'Residuals for the quadratic fit'

Higher order polynomials - are more terms always better?

Polynomial modeling with polyfit is indeed simple and easy to do. In fact, we might decide to just use higher and higher order polynomials, always chasing the term we truncated in the previous model. After all, if a quadratic model is better than a linear one, then why not go to a cubic? Ten terms must be better than nine. Look at a tenth order model.

P10 = polyfit(x,y,10);
disp(P10(1:6))
disp(P10(7:11))

    0.0000    0.0000    0.0000    0.0002    0.0014    0.0083
    0.0417    0.1667    0.5000    1.0000    1.0000

We can write the Maclaurin series representation for the exponential function as

We can compare P10 to the coefficients of the known series. How well did we do in the fit?

series = 1./factorial(10:-1:0);
disp(series(1:6))
disp(series(7:11))

    0.0000    0.0000    0.0000    0.0002    0.0014    0.0083
    0.0417    0.1667    0.5000    1.0000    1.0000

Note that the higher order coefficients deviate somewhat from the known series, although the lower order terms appear to be quite accurate.

yhat = polyval(P10,x);
plot(x,y,'bo',x,yhat,'r-')
xlabel X
ylabel Y
grid on
title 'Tenth order polynomial fit'

The residuals oscillate tightly around zero. In fact, they are so small that this last polynomial begins to approach a true interpolating polynomial. Perhaps if we just add a few more terms we may get there. The numerical issues of floating point arithmetic will often preclude true interpolation down to the least significant bit anyway.

res = y - yhat;
plot(x,res,'bo')
xlabel X
ylabel Residuals
grid on
title 'Residuals for the tenth order fit'

What do you do with interpolation?

I'll start talking about true interpolation in my next blog. But remember that interpolation is different from the approximations provided by polyfit or any other regression modeling tool.

Until that time, please give me your comments on this blog, or ideas for future blog topics on interpolation or modeling in general. Do you have some interesting applications of interpolation? Some interesting problems?

Get the MATLAB code

Published with MATLAB® 7.6