INTRODUCTION TO MATLAB
Chapter 11
Interpolation and Extrapolation
Since Matlab only represents functions as arrays of values a common problem that
comes
up is finding function values at points not in the arrays. Finding function values
between
data points in the array is called interpolation, finding function values beyond
the endpoints
of the array is called extrapolation. A common way to do both is to use nearby
function
values to define a polynomial approximation to the function that is pretty good
over a small
region. Both linear and quadratic function approximations will be discussed
here.
11.1 Linear Interpolation and Extrapolation
A linear approximation can be obtained with just two data points, say (x_{1}, y_{1})
and (x_{2}, y_{2}).
You learned a long time ago that two points define a line and that the two-point
formula
for a line is
This formula can be used between any two data points to linearly interpolate.
For example,
if x in this formula is half-way between x_{1} and x_{2} at x = (x_{1}+x_{2})/2 then it is
easy to show
that linear interpolation gives the obvious result y = (y_{1} + y_{2})/2.
But you must be careful when using this method that your points are close enough
together to give good values. In Fig. 11.1, for instance, the linear
approximation to the
curved function represented by the dashed line "a" is pretty poor because the
points x = 0
and x = 1 on which this line is based are just too far apart. Adding a point in
between at
x = 0.5 gets us the two-segment approximation "c" which is quite a bit better.
Notice also
that line "b" is a pretty good approximation because the function doesn't curve
much.
This linear formula can also be used to extrapolate. A common way extrapolation
is
often used is to find just one more function value beyond the end of a set of
function pairs
equally spaced in x. If the last two function values in the array are f_{N-1} and f_{N}, it is easy
to show that the formula above gives the simple rule
which was used in the Matlab code in Sec. 10.1
Figure 11.1 Linear interpolation only works well over intervals
where the function is straight.
You must be careful here as well: segment "d" in Fig. 11.1 is the linear
extrapolation of
segment "b", but because the function starts to curve again "d" is a lousy
approximation
unless x is quite close to x = 2.
11.2 Quadratic Interpolation and Extrapolation
Quadratic interpolation and extrapolation are more accurate than linear because
the
quadratic polynomial ax^{2} + bx + c can more easily t curved functions than the
linear
polynomial ax + b. Consider Fig. 11.2. It shows two quadratic fits to the curved
function.
The one marked "a" just uses the points x = 0, 1, 2 and is not very accurate
because these
points are too far apart. But the approximation using x = 0, 0.5, 1, marked "b",
is really
quite good, much better than a two-segment linear t using the same three points
would
be.
Unfortunately, the formulas for quadratic fitting are more difficult to derive. (The
Lagrange interpolation formulas, which you can find in most elementary numerical
analysis
books, give these formulas.) But for equally spaced data in x, Taylor's theorem,
coupled with
the approximations to the first and second derivatives discussed in the section on
numerical
derivatives, make it easy to derive and use quadratic interpolation and
extrapolation. We
want to do it this way because it uses the approximate derivative formulas we
used in
Sec. and illustrates a technique which is widely used in numerical
analysis.
You may recall Taylor's theorem that an approximation to the function f(x) near
the
point x = a is given by
Let's use this approximation (ignoring all terms beyond the quadratic term in (x
- a))
near a point (x_{n}, f_{n} ) in an array of equally spaced x values. The grid spacing
in x is h.
Figure 11.2 Quadratic interpolation follows the curves better if
the curvature doesn't change sign.
An approximation to Taylor's theorem that uses numerical derivatives in this
array is then
given by
This formula is very useful for getting function values that aren't in the
array. For instance,
it is easy to use this formula to obtain the interpolation approximation to f(x_{n}
+ h/2)
and also to find the quadratic extrapolation rule
11.3 Interpolating With polyfit and polyval
You can also use Matlab's polynomial commands (which you learned about in
Chapter 8)
to build an interpolating polynomial. Here is an example of how to use them to
find a
5th-order polynomial t to a crude representation of the sine function.
Example 11.3a (ch11ex3a.m) % Example 11.3a (Physics 330) clear, close all, dx=pi/5, plot(x,y,'b*',xi,yi,'r-',xi,sin(xi),'c-') % display the difference between the polynomial fit and |
11.4 Matlab Commands Interp1 and Interp2
Matlab has its own interpolation routine interp1 which does the things discussed
in this
section automatically. Suppose you have a set of data points {x, y} and you have
a different
set of x-values {x_{i}} for which you want to find the corresponding {y_{i}} values by
interpolating
in the {x, y} data set. You simply use any one of these three forms of the
interp1 command:
yi=interp1(x,y,xi,'linear')
yi=interp1(x,y,xi,'cubic')
yi=interp1(x,y,xi,'spline')
Here is an example of how each of these three types of interpolation works on a
crude
data set representing the sine function.
Example 11.4a (ch11ex4a.m)
% Example 11.4a (Physics 330) % make the crude data set with dx too big for % interpolate on the coarse grid to % linear interpolation plot(x,y,'b*',xi,yi,'r-') % cubic interpolation % plot the data and the interpolation % spline interpolation % plot the data and the interpolation |
Matlab also knows how to do 2-dimensional interpolation on a data set of
{x, y, z} to
find approximate values of z(x, y) at points {x_{i}, y_{i}} which don't lie on the data
points {x, y}.
You could use any of the following
zi = interp2(x,y,z,xi,yi,'linear')
zi = interp2(x,y,z,xi,yi,'cubic')
zi = interp2(x,y,z,xi,yi,'spline')
This will work fine for 1-dimensional data pairs {x_{i}, y_{i}}, but you might want to
do this
interpolation for a whole bunch of points over a 2-d plane, then make a surface
plot of
the interpolated function z(x, y). Here's some code to do this and compare these
three
interpolation methods (linear, cubic, and spline).
Example 11.4b (ch11ex4b.m)
% Example 11.4b (Physics 330) % build the full 2-d grid for the crude x and y data [X,Y]=meshgrid(x,y), %********************************************************* % Note that because the interpolation is linear the mesh is finer xf=-3:.1:3, %********************************************************* ZF=interp2(X,Y,Z,XF,YF,'cubic'), %********************************************************* ZF=interp2(X,Y,Z,XF,YF,'spline'), |
For more detail see Mastering Matlab 6, Chapter 19.