Philip Gage is a software engineer in Colorado Springs, currently developing Ada graphics software with X Windows. He has a BS degree in computer science from the University of Colorado and eleven years programming experience developing FORTRAN, BASIC, Pascal, C, and Ada programs on a variety of platforms. His programming interests include data compression, algorithms, graphics, and compiler design. He has written an article, "A New Algorithm for Data Compression," in the February 1994 issue of The C Users Journal.
Introduction
Interpolation is the process of estimating values between known values, and has many practical uses. The same process is called extrapolation if the estimated value lies outside the range of known values.This article describes an algorithm that generates interpolation functions using Newton's method of divided differences. The algorithm inputs a set of (x,y) data points and outputs a polynomial function y = f(x) of any desired degree which defines a curve through or near the points.
The method of divided differences can be extended to generate other useful classes of functions, including the exponential and factorial. These techniques form a powerful set of tools for interpolation, extrapolation, curve fitting, and numerical analysis. Such techniques are well known to mathematicians but may be new to most programmers.
The method of finite differences was first described centuries ago by the famous mathematician Sir Isaac Newton. The Difference Engine, an early computer built by Charles Babbage and programmed by Lady Ada Lovelace, used the difference method to print tables of values for many kinds of functions.
Differences can also be used for differentiation, integration, correcting errors in tables, and other purposes. Such difference techniques were widely used for their power and efficiency before modern computers became available.
Most of us have faced this type of problem, "Given the numbers 1, 3, 6, 10, what is the next number in the series ?" After some thought, we realize that the differences between the numbers are increasing 2, 3, 4. Therefore, the next difference should be 5 and the next number should be 15.
This example illustrates a profound mathematical law, that subtracting adjacent values can reduce data complexity to a simple relationship that allows backtracking to fill in additional data values. This is actually a form of numerical differentiation, but you don't need to know any calculus to understand the method, just a little algebra.
Polynomials
The method of differences is associated with an important class of functions called polynomials. Recall from high school algebra that a polynomial y = f(x) has the following form:
c0 + c1x + c2x2 + c3x3+ ... + cnxnwhere c0, c1, ... are constant numbers called coefficients. The degree of each term is the exponent of x, and the degree of the polynomial is that of the highest degree term.The first difference (usually denoted by D, the Greek letter delta) of a set of (x,y) values with evenly spaced x coordinates is obtained by subtracting adjacent y values. Using C array notation, the y values are y[0], y[1], ...
Higher-order differences are calculated similarly by subtracting previous differences, as shown for the function y = x2 in Table 1. Each entry in a triangular matrix is computed as the difference of the two entries to its left, starting with the original y values as the zeroth difference.
First difference of y[n]:
y[n+1] - y[n]Second difference of y[n]:
(y[n+2] - y[n+1]) - (y[n+1] - y[n])Assuming the set of data points follows a polynomial distribution, each difference operation reduces the degree of the polynomial by one. The first difference of an n-th degree polynomial is a polynomial of degree n--1, the second difference is a polynomial of degree n--2, and so on. The n-th difference is a constant, and all higher differences are zero. The difference interpolation method thus has the advantage that an underlying polynomial distribution in the data is easily detected by the presence of a constant difference.
Babbage's Method
Additional (x,y) values can be derived easily from a difference table by reversing the process. In Table 1, the second difference is a constant 2. Therefore, the next first difference value is 9 + 2 = 11. The next y value is 25 + 11 = 36. This is precisely the method used by Babbage's Difference Engine to print tables of polynomial values efficiently, using only addition.If the n-th difference of a data set is constant, the data can be represented by a polynomial of degree n. Every set of n+1 points can be generated by a polynomial with degree of at most n, since the last difference column has only one entry and is therefore a constant difference.
The n-th degree polynomial term can be determined from this constant difference. All other terms of degree n--1 and lower have been reduced to zero and filtered out of the data. Only the effect of the highest order term remains. If the n-th difference is a constant c and the x values are uniformly spaced a distance h apart, the high-order term is given by:
(cxn)/(hnn!)Once the high-order term is identified, it can be removed from the data set by evaluating the term at each x point and subtracting it from the corresponding y value, reducing the polynomial order of the data by one. This process is then repeated, solving each high order term and removing it from the data, until the y values are reduced to zero and the complete solution is known. There are other ways of determining the interpolating polynomial from a table of differences, but this method is simple and intuitive.To illustrate the overall solution method, consider the arithmetic series 1, 3, 6, 10. Assuming x values beginning at 1, we can build a simple difference table. The first differences are 2, 3, 4, ... Thus, the second difference is a constant 1, and the high-order term is given by:
(cxn) / (hnn!) = (1x2) / (122!) = 0.5x2Evaluating the term 0.5 x2 at each x value and subtracting from the original y values gives a difference table with a constant first difference of 0.5. So the high-order term is given by:
(cxn) / (hnn!) = (0.5x1) / (111!) = 0.5xSubtracting 0.5 x from the remaining y values leaves all zeros, so the complete solution is
0.5x2 + 0.5xDivided Differences
One limitation of the difference method is that the x values must be ordered and evenly spaced. This restriction can be removed by using divided differences. If each difference value is divided by the x interval it spans, the result is called a divided difference. See Table 2. First divided difference of y[n]:
(y[n+1] - y[n]) / (x[n+1] - x[n])Second divided difference of y[n]:
((y[n+2] - y[n+1]) - (y[n+1] - y[n])) / (x[n+2] - x[n])If the n-th divided difference is a constant c, the high-order term is c xn.
Exponentials
The method of differences can be extended from polynomial to exponential functions by a process I call the method of quotients, using division instead of subtraction. Similar to polynomials, but using different operators, this class of exponential functions has the form:
c0 * c1x * c2x2 * c3x3 * ... * cnxnThe first quotient is obtained by dividing adjacent y values. Higher-order quotients are calculated similarly by dividing previous quotients. Once a constant value is reached, all further quotients will be 1, as illustrated in Table 3.First quotient of y[n]:
y[n+1] / y[n]Second quotient of y[n]:
(y[n+2] / y[n+1]) / (y[n+1] / y[n])If the n-th quotient is a constant c, with x values evenly spaced a distance h apart, the high order term is given by
c(xn)/(n!hn)Each high-order term is removed from the y data values by division, instead of subtraction as with polynomials. When the data values are reduced to 1, the complete solution is known. The quotient method has the disadvantage that y values must be either all positive or all negative, since it is based on exponentials. If necessary, a constant can be manually added to the y values to normalize them to a positive or negative range before applying the algorithm.In the same manner as divided differences, a method I call root quotients can be used to remove the limitation that the x coordinates be uniformly spaced. Instead of division, the operation used is taking a root, as shown in Table 4. Taking the n-th root is equivalent to raising to the 1/n power.
First root quotient of y[n]:
(y[n+1] / y[n]) ^ 1/(x[n+1] - x[n])If the n-th root quotient is a constant c, the high order term is given by
cxnOther functions
Besides polynomial and exponential, other classes of functions can be handled by similar methods. Performing one quotient operation followed by differences solves a class of factorial expressions. Doing differences before quotients handles the sum of an exponential and a polynomial function, by filtering out the polynomial terms without affecting the exponential terms.These methods can also be extended to handle functions of more than one variable, such as z = f(x,y), by a process equivalent to partial differentiation. This requires using multi-dimensional arrays for each difference or quotient, but is straightforward to program.
Implementation
Pseudocode for the divided-difference and root quotient algorithms is shown in Listing 1 and Listing 2. Listing 3 contains an interactive C program which generates both polynomial and exponential interpolation functions using these algorithms. For simplicity, error checking is minimal. You may want to add checks for dividing by zero and pow function errors. Listing 4 shows the output from a typical run.The program prompts the user for a set of (x,y) points, asks whether a polynomial or exponential function is wanted, and requests the desired degree of the solution. Output is the interpolation function and a table of results comparing the original and interpolated values.
To generate an exact solution which goes through all n points, use degree n--1. Lower-degree solutions will be approximations with some error, but may be accurate enough for many applications, such as curve fitting. To improve low-degree approximations, each coefficient is computed using the average of all the difference values assumed to be constant.
The maximum number of data points is set to only ten since high-order interpolation functions tend to oscillate wildly and can produce bad results. Large data sets should be split into smaller ranges and interpolated separately. However, the difference and quotient algorithms cannot join segments together smoothly like some spline techniques.
Function interpolate implements the divided-difference and root-quotient algorithms and returns an array of coefficients for the generated function. Since the code for both methods is very similar, it is combined, using a polynomial/exponential flag to control execution of the few expressions that differ. This program does not need to check for the appearance of a constant set of difference values, but simply generates a polynomial of the requested degree (or lower, since some coefficients may be zero).
Function evaluate uses Horner's rule to evaluate a polynomial or exponential function at any x value, using the coefficient array returned by interpolate. This function can be called by applications to interpolate new data points.
There are two print functions to print results for informational purposes, print_function prints the symbolic interpolation function term by term, while print_table prints a table of results comparing the original data with the interpolated values.
Summary
The divided-difference algorithm is an effective method for generating polynomial solutions or approximations to a set of data points. The root quotient method works similarly for exponential functions. These methods can be extended to other classes of functions, such as factorial and functions of more than one variable, providing a flexible methodology for interpolation and extrapolation.
References
1. Levy, H., Lessman, F., Finite Difference Equations, Macmillan 1961.2. Maron, M., Numerical Analysis, Macmillan 1982.