Interpolation is a mathematician's version of "connect the dots." More formally, it is a way of approximating the value of a function at various points using a limited number of values known at others. The ability to do this is important in many applications. For example, imagine a program that requires the use of wind velocities predicted for various locations and altitudes around the world. Here it would be prohibitively costly (and wasteful) to store values for every location and altitude required. Instead, you just need a good idea of the winds at a limited number of points; you can then interpolate between them. Another use of interpolation is to replace more complicated functions with simpler ones so that operations such as differentiation and integration become easier. Human beings are perhaps the best interpolators of all. If you flash a sequence of static images before your eyes fast enough (say, of a horse running), your brain will fill in the gaps.
This article presents five methods for interpolating functions using a computer. Each method is based on the construction of polynomials. Polynomials are an important part of the interpolation process because they have a simple structure, which makes them easy to work with in a variety of contexts. Specifically, this article presents interpolation using Lagrange polynomials, interpolation using Newton polynomials, piecewise-linear interpolation, piecewise-cubic Hermite interpolation, and piecewise-cubic Bessel interpolation. The article concludes with a brief discussion of techniques for estimating interpolation error.
Interpolation Using Lagrange Polynomials
The goal of polynomial interpolation is as follows: given a set of n + 1 points x0,...,xn on the interval I = [a, b] over which we know the value of a function f(x), construct a polynomial that satisfies the equation below.
The polynomial that satisfies this is said to interpolate the function f(x). For a set of n + 1 points, this polynomial is of degree <= n. The importance of p(x) is that once you have found it, you can use it to approximate f(x) at values of x not in x0,...,xn, provided these values are on the interval I. The standard way to represent a polynomial in mathematical discussions is in the power form, shown below. This is generally the most familiar form. Provided an is nonzero, this polynomial is said to be of degree n.
Although the power form of a polynomial is common, other polynomial forms are more useful in some contexts. For example, the Lagrange form, shown below, is particularly useful in the context of interpolation.
where
For a set of n + 1 points x0,...,xn along a function f(x), you would construct the Lagrange form of a polynomial with n = 2 as follows:
Since each function lk(x) is a product of n linear factors, each lk(x) is a polynomial of exact degree n; consequently, the Lagrange form as a whole describes a polynomial of degree n.
Considering the Lagrange form further, notice that lk(x) becomes 0 at xi for all i
k, and 1 when i = k. Therefore, you can make the statement below.
In other words, the coefficients a0,...,an in the Lagrange form are the values of the polynomial p(x) at the points x0,...,xn. Replacing a0,...,an with the values of f(x) at x0,...,xn yields the following formula, which produces a polynomial of degree <= n that interpolates f(x) at x0,...,xn.
This is the Lagrange formula for the interpolating polynomial. Listing 1 shows how to use this formula in a program.
Example 1
This example uses Lagrange polynomials to interpolate a function f(x) describing the lift produced by the wing of a large commercial aircraft as a function of airspeed at a fixed altitude and wing configuration. Suppose you know the values of f(x) at x0 = 100 ft/sec, x1= 500 ft/sec, and x2 = 900 ft/sec, which are as follows:
f(100) = 7882.1 lbs f(500) = 197052.5 lbs f(900) = 638450.1 lbs
To compute the lift produced at x = 600 ft/sec, for example, you perform the calculations shown below. Compute the lift for other values of x in a similar manner.
Figure 1 shows the results of approximating f(x) at x = 200, x = 400, x = 600, and x = 800 ft/sec using the three points known for f(x) presented previously.
In this example, since there are three points at which f(x) is known, the interpolating polynomial is of degree <= 2. It turns out that the function for lift is, in fact, a second-degree polynomial when written in terms of airspeed. Therefore, Figure 1 is a good representation of f(x). The discussion on error later presents a more quantitative explanation of why this is the case.
Interpolation using Lagrange polynomials is a nice method in such cases when you know exactly how many points you need to use to get a good impression of f(x). However, in practice, you often don't know this. In these situations, it is common to compute several interpolating polynomials pk(x) of increasing degree <= k (i.e., p0(x), p1(x), p2(x), etc.) until |pk(x) - pk-1(x)|(k >= 1) is less than some tolerance. When interpolating using Lagrange polynomials, you must compute each successive pi(x) from scratch, which ends up being more work than necessary. Interpolation using Newton polynomials offers a more efficient approach.
Interpolation Using Newton Polynomials
In the previous section, you saw how the Lagrange form is one alternative to the power form for representing a polynomial. Another polynomial form relevant to interpolation is the Newton form, shown below. Provided an is nonzero, this polynomial is said to be of degree n.
The variables x0,...,xn-1 are centers, and in the case of polynomial interpolation, are the first n of the n + 1 interpolation points. To construct the interpolating polynomial in the Newton form, replace each ak with f[x0,...,xk], k = 0,...,n.
This is the Newton formula for the interpolating polynomial, which is often written more succinctly in the form that follows. Each f[x0,...,xk] is the kth divided difference of f(x) at points x0,...,xk.
To compute the value of each f[x0,...,xk], the coefficients of the interpolating polynomial, use the procedure described next. Let f[x0] = f(x0) and f[x0, x1] = (f(x1) - f(x0)) / (x1 - x0), which you can verify algebraically by manipulating p1(x) derived from the Lagrange formula as outlined below.
Compute higher-order divided differences using the formula below, which assumes n + 1 interpolation points. Notice how each divided difference depends on computing others before it.
The best way to perform these calculations in an orderly fashion is with a divided-difference table. A divided-difference table consists of several rows. In the top row, place the values x0,...,xn. In the second row, place the values f[x0] = f(x0),...,f[xn] = f(xn). Then, to compute each divided difference in the remainder of the table, draw a diagonal from each divided difference back to f[xi] and f[xi+k]. To get xi and xi+k in the denominator, proceed straight up from f[xi] and f[xi+k]. The two divided differences in the numerator are the two immediately above the one being computed. Once you have populated the entire table, the coefficients for the interpolating polynomial are the divided differences at the far left of each row, in order, beginning with the second row.
Figure 2 illustrates a divided-difference table generated for four points. As mentioned earlier, one of the nice things about interpolation using Newton polynomials is that it lets you proceed more efficiently than the Lagrange form when computing polynomials pk(x) of increasing degree <= k until |(pk(x) - pk-1(x)|(k >= 1) is less than some tolerance. To understand this, let qk(x) be the sum of the first k + 1 terms of pn(x) for some k between 0 and n. In this case, each of the remaining terms has the factor (x - x0)...(x - xk); thus, you can write pn(x) in the form that follows, where r(x) is a polynomial of no interest for now.
Since (x - x0)...(x - xk)r(x) becomes 0 when x is any point in x0,...,xk, qk(x) already interpolates f(x) at x0,...,xk. Since qk(x) is a polynomial of degree <= k, it follows that qk(x) = pk(x); therefore, qk(x) is the polynomial of degree <= k that interpolates f(x) at x0,...,xk. This is important because it means that you can construct each successive pk(x) incrementally from its predecessor, pk-1(x). To compute pk(x), add pk-1(x) to the next term in the Newton form, as shown below.
Listing 2 shows how to perform this process in a program. Note that when compared with interpolation using Lagrange polynomials, interpolation using Newton polynomials usually requires computing a polynomial one degree larger than what you might expect; consequently, one additional interpolation point is required as well. The polynomial of one degree larger is pk(x), which must be used to determine when |pk(x) - pk-1(x)| is less than the tolerance.
Example 2
This example uses Newton polynomials to approximate values of f(x) = sin x given nine interpolation points: xi = ip/4, i = 0,...,8, as follows:
f(0) = 0.0 f(p/4) = 0.7071068 f(p/2) = 1.0 f(3p/4) = 0.7071068 f(p) = 0.0 f(5p/4) = -0.7071068 f(3p/2) = -1.0 f(7p/4) = -0.7071068 f(2p) = 0.0
To compute the value of f(x) at x = 3p/8 = 1.17810 so that |pk(x) - pk-1)(x)| < 0.001 (k >= 1), for example, you end up computing p0(x),...,p6(x). This process is illustrated below.
...
Interpolation using Newton polynomials is nice in cases when you know that f(x) is simple enough that it can be approximated to within the desired tolerance using a relatively small number of interpolation points (say <= 10). However, since interpolation using Newton polynomials (as well as Lagrange polynomials) involves constructing a single polynomial to interpolate f(x) at each of the interpolation points across the entire interval of interest, increasing the number of points complicates the overall interpolation process. When interpolating more complicated functions, such as functions involving frequent changes in concavity, more points are required. In these cases, one of the following piecewise approaches is preferable.
Piecewise-Linear Interpolation
To interpolate a function in a piecewise manner, you break it logically into a series of pieces and interpolate each piece individually. In piecewise-linear interpolation, you interpolate each piece using a first-degree polynomial. This means that two points are required to interpolate each piece, and the curve between each pair of points is a straight line (hence the term "linear"). Formally, given a series of n + 1 values for f(xi), i = 0,...,n, x0 < ... < xn, you approximate f(x) at some point x by determining the interval [xk, xk+1] on which x resides and performing the calculation below using xk and xk+1 as the interpolation points.
The advantage of piecewise-linear interpolation is that although f(x) may be complicated, the polynomials that interpolate the pieces of f(x) are not. Furthermore, once you understand interpolation using Newton polynomials, it's easy to see that piecewise-linear interpolation is just a simplified case, applied several times instead of just once.
Listing 3 shows how to perform piecewise-linear interpolation in a program. As an illustration, Figure 3a depicts the first second of ground motion associated with a minor earthquake recorded by the Berkeley Digital Seismic Network. The tremor occurred in Woodside, CA (near San Francisco) as a result of one of the largest eucalyptus trees in the United States toppling to the ground. Figure 3b shows the linear pieces that result from interpolating the curve in Figure 3a using a linear-form interpolating polynomial with interpolation points equally spaced at xi = 0.05i seconds, i = 0,...,20. One thing you may notice in Figure 3b is that as the distance between interpolation points decreases, the lines that interpolate the pieces of f(x) decrease as well; consequently, the error associated with each piece becomes smaller. The discussion on error later presents a more quantitative explanation of why this is the case. Another way to reduce the error is to use higher-order polynomials to interpolate the pieces. Cubic polynomials are particularly popular.
Piecewise-Cubic Hermite Interpolation
Piecewise-cubic Hermite interpolation is similar to piecewise-linear interpolation. It too involves breaking a function logically into a series of pieces that are then interpolated individually. However, in piecewise-cubic Hermite interpolation, as in all piecewise-cubic methods, you interpolate each piece using a third-degree polynomial. Four points are thus required to interpolate each piece; consequently, the resulting curve fits the piece closer than a straight line would in most cases.
Formally, to approximate f(x) at some point x, given a series of n + 1 values for f(xi), i = 0,...,n, x0 < ... < xn, you determine the interval [xk, xk+1] on which x resides and evaluate the third-degree interpolating polynomial. Considering this, if xk and xk+1 are the only two points for which you know values for a piece of f(x), it's worth asking where the other two of the four requisite points come from. Different piecewise-cubic interpolation methods answer this question in different ways. Piecewise-cubic Hermite interpolation answers the question by selecting the points xk, xk, xk+1, and xk+1. Selecting these points ensures two things: first, that p(x) agrees with f(x) at x = xk and x = xk+1 (as with any interpolating polynomial); and second, that p'(x), the derivative of p(x), also agrees with f'(x), the derivative of f(x), at x = xk and x = xk+1. The Newton form of the interpolating polynomial generated from these points is shown below.
In this polynomial, each of the four coefficients, call them a0, a1, a2, and a3, is a divided difference expressed in terms of xk and xk+1. To determine the value of each coefficient, use the formulas that follow. To see where these formulas come from, construct a divided-difference table using the points x0 = xk, x1 = xk, x2 = xk+1, and x3 = xk+1. In the table, wherever you see (f(xk+1) - f(xk+1))/(xk+1 - xk+1), replace it with f'(xk+1); wherever you see (f(xk) - f(xk))/(xk - xk), replace it with f'(xk). Listing 4 shows how to perform piecewise-cubic Hermite interpolation in a program.
Notice that piecewise-cubic Hermite interpolation requires the value of f'(x) at each interpolation point in addition to the value of f(x). In practice, this requirement for f'(x) is sometimes problematic. You may be able to get values for f'(x) when interpolating a function whose values have been stored in a function table; however, there are many cases when obtaining these values will be either difficult or not possible at all. When this is true, you can approximate f'(x). This is called piecewise-cubic Bessel interpolation.
Piecewise-Cubic Bessel Interpolation
Piecewise-cubic Bessel interpolation is similar to piecewise-cubic Hermite interpolation, except that it substitutes an approximation for f'(x) in the computation of coefficients for the interpolating polynomial. In piecewise-cubic Bessel interpolation, you carry out the same computations as those presented for a0, a1, a2, and a3 above, except you use the formula below for each f'(xk). This formula requires the use of two additional points, one to the left of x0 and one to the right of xn, to give some number for derivatives when interpolating the leftmost and rightmost pieces. Listing 4 performs piecewise-cubic Bessel interpolation when you use approximations for f'(x).
Estimating Interpolation Error
This section describes the error associated with the various interpolation methods presented in this article. Fundamentally, the interpolation error en(x) associated with any interpolating polynomial pn(x) of degree <= n interpolating a function f(x) is defined by the equation below.
Of course, normally you won't know f(x), so you must look for a way to estimate the value of en(x) in practice. Assuming X to be some point different from x0,...,xn, let pn+1(x) be the polynomial that interpolates f(x) at x0,...,xn and X as well. In other words, pn+1(X) = f(X). Recall that you can construct the interpolating polynomial pn+1(x) from the interpolating polynomial of one degree smaller, pn(x), using the Newton form. Therefore, since pn+1(x) interpolates X, the error in pn(x) is given by the next term added to pn(x), which produces pn+1(x).
This still doesn't help you much, though, because to compute f[x0,...,xn, X], the need for f(x) remains. However, it turns out that f[x0,...,xn, X] is closely related to the (n + 1)st derivative of f(x). Therefore, when you know a bound for f(n + 1)(x), call it M, over the entire interval you are interpolating, you can use it to obtain a bound (albeit a broad one) for the error. Considering this, the formula below is used to estimate interpolation error.
The variable x in the formula is a rarely-known value that depends on the point X at which you desire the error estimation. It appears in the formula because interpolation error varies depending on where you look along f(x). This makes sense when you think of how an interpolating polynomial often yields approximations for f(x) that are better at some points than others (which is why I mention Chebyshev points below). Since you almost never know x, the key is to replace f(n+1)(x) with the bound M mentioned a moment ago.
One interesting thing to notice in the error formula is that when you use n + 1 points to interpolate a function (thus producing an interpolating polynomial of degree n), f(n+1)(x) becomes 0; consequently, the interpolation error becomes 0 as well, regardless of the value of x. This explains quantitatively why when interpolating the second-degree function for lift in the first example, for instance, interpolation with three points produced a good representation of the second-degree function.
From the basic error formula, you can derive formulas for piecewise-linear interpolation and piecewise-cubic Hermite interpolation as well, which are listed in order below (assuming interpolation over the interval I = [a,b] and n + 1 interpolation points). There are two things to notice in these formulas: first, you can reduce the interpolation error as small as necessary by making the distance between interpolation points smaller for all k; and second, you can reduce the error by using interpolating polynomials of a greater degree for each of the pieces (as in the case of trying piecewise-cubic interpolation as opposed to piecewise-linear interpolation).
Another aspect of interpolation error is how uniformly well a polynomial interpolates a function across an entire interval of interest. This depends on the placement of the interpolation points. For this, I encourage you to research the significance of Chebyshev points. The reference [1] provides an excellent description.
Summary
This article has presented several types of interpolation methods along with brief explanations on when to apply them. To review, interpolation using Lagrange polynomials is a nice way to understand the existence of the interpolating polynomial, and is useful when you know the number of points required to give a good impression of a function. Interpolation using Newton polynomials is a good approach when you need to compute polynomials of increasingly larger degree to find one that approximates a function within a desired tolerance. Piecewise approaches (piecewise-linear interpolation, piecewise-cubic Hermite interpolation, and piecewise-cubic Bessel interpolation) are good for interpolating more complicated functions (functions that require the use of more than ten interpolation points). All piecewise approaches work by logically breaking functions into pieces that are interpolated individually. Piecewise approaches vary in the degree of the polynomials used to interpolate the pieces as well as in which points are selected as interpolation points.
You can use the interpolation methods presented in this article to solve interpolation problems found in many applications. For a more extensive treatment, there is a great deal of research and academic literature on the subject to which you can refer. [1], [2], and [3] are some good examples.
References
[1] Samuel D. Conte and Carl de Boor. Elementary Numerical Analysis: An Algorithmic Approach (McGraw-Hill, 1980). This book provides a particularly good mathematical presentation of interpolation methods.
[2] William H. Press, Saul A. Teukolsky, William T. Vetterling, and Brian P. Flannery. Numerical Recipes in C: The Art of Scientific Computing (Cambridge University Press, 1993). This is a popular reference regarded by many engineers, scientists, and software developers as particularly friendly. It also has the benefit of being available in its entirety online at http://www.ulib.org/webRoot/Books/Numerical_Recipes.
[3] Michael C. Kohn. Practical Numerical Methods: Algorithms and Programs (Macmillan Publishing Co., 1987). This book provides nice illustrations for visualizing how some of the methods presented in this article work.
Kyle Loudon currently works for Pixo, Inc. in Cupertino, CA, a company developing web technologies, application software, and a graphics-rich user interface for mass-market wireless phones. Before Pixo, Kyle worked as the technical leader of graphical user interface development at Jeppesen Dataplan and as a system programmer for IBM. He has a B.S. in computer science from Purdue University and has done advanced study at Stanford University. He is the author of Mastering Algorithms with C, published by O'Reilly & Associates.