A New Method of Interpolation and S m o o t h Curve Fitting Based on Local Procedures

HIROSHI A K I M A

ESSA Research Laboratories,* Boulder, Colorado ABSTRACT. A new mathematical method is developed for interpolation from a given set of data points in a plane and for fitting a smooth curve to the points. This method is devised in such a way that the resultant curve will pass through the given points and will appear smooth and natural. It is based on a piecewise function composed of a set of polynomials, each of degree three, at most, and applicable to successive intervals of the given points. In this method, the slope of the curve is determined at each given point locally, and each polynomial representing a portion of the curve between a pair of given points is determined by the coordinates of and the slopes at the points. Comparison indicates that the curve obtained by this new method is closer to a manually drawn curve than those drawn by other mathematical methods. KEYWORDSAND PHRASES: interpolation, polynomial, slope of curve, smooth curve fitting
CR CATEGORIES:

5.13

1. Introduction To determine a relation b e t w e e n two variables, we either perform c o m p u t a t i o n s or make measurements. T h e result is given as a set of discrete d a t a points in a plane. Knowing t h a t the relation can be represented b y a s m o o t h curve, we next t r y to fit a smooth curve to the set of d a t a points so t h a t it will pass t h r o u g h all the points. Manual drawing is t h e m o s t primitive m e t h o d for this purpose and results in a reasonable curve if it is done b y a well-trained scientist or engineer. B u t , since it is very tedious and time consuming, we wish to let a c o m p u t e r draw t h e curve. T h e computer m u s t then be provided with necessary instructions for m a t h e m a t i c a l l y interpolating additional points between the given d a t a points. There are several m a t h e m a t i c a l m e t h o d s of interpolating a single-valued function from a given set of values [3, 4, 6], b u t their application to curve fitting sometimes results in a curve t h a t is v e r y different f r o m one d r a w n manually. T h e c o m m o n difficulty is t h a t the resultant curve sometimes shows u n n a t u r a l wiggles. This seems inevitable if we m a k e a n y assumption concerning t h e functional form for t h e whole set of given d a t a points o t h e r t h a n the continuity and t h e smoothness of the curve. W h e n we t r y to fit a s m o o t h curve manually, we do n o t assume a n y functional form for t h e whole curve; we draw a portion of t h e curve based on a relatively small number of points, w i t h o u t t a k i n g into account t h e whole set of points. This local aspect is a v e r y i m p o r t a n t feature of m a n u a l curve fitting and is t h e basis for our new method. I t is n o t easy to develop a m a t h e m a t i c a l m e t h o d of s m o o t h curve fitting based on the local procedure. If one of t h e existing m a t h e m a t i c a l m e t h o d s is applied locally or piecewise w i t h o u t special consideration, t h e continuity of the function or * Institute for Telecommunication Sciences

Journal of the Association for Computing Machinery, Vol. 17, No. 4, October 1970, pp. 589-602.

590

HIROSHI

AKIMA

its first-order derivative at the junction points cannot be generally guaranteed. However, if we can determine the slope of the curve at each given point locally, we obtain a smooth curve by piecewise application of one of the existing mathematical methods. The osculatory interpolation method [1, 5], which is based on a piecewise function composed of a set of third-degree polynomials with slopes at the junction points locally determined, was developed along this line. But, sometimes this method also gives unnatural curves. Recently the author proposed a method of interpolation and smooth curve fitting t h a t is based on a piecewise function with slopes at the junction points locally determined under a geometrical condition [2]. In the present study this method is further developed with an improved condition for determining the slope of the
curve.

In this paper, except in Appendix B, we describe a new method of interpolation and smooth curve fitting that is applicable to a single-valued function. However, the basic idea of our new method is also applicable to a multiple-valued function, and a method applicable to this function is outlined in Appendix B.
2. N e w Method

Our method is based on a piecewise function composed of a set of polynomials, each of degree three, at most, and applicable to successive intervals of the given points. We assume that the slope of the curve at each given point is determined locally by the coordinates of five points, with the point in question as a center point, and two points on each side of it. This is discussed in the next section. A polynomial of degree three representing a portion of the curve between a pair of given points is determined by the coordinates of and the slopes at the two points. This interpolation procedure is described in Section 2.2. Since the slope of the curve must thus be determined also at the end points of the curve, estimation of two more points is necessary at each end point. This estimation procedure is described in Section 2.3. 2.1. SLOPE OF THE CURVE. With five data points 1, 2, 3, 4, and 5 given in a plane, we seek a reasonable condition for determining the slope of the curve at point 3. It seems appropriate to assume t h a t the slope of the curve at point 3 should approach t h a t of line segment 2-3 when the slope of ~ approaches t h a t of ~ . It is also highly desirable that the condition be invariant under a linear-scale transformation of the coordinate system. With these rather intuitive reasonings as a guideline, the condition of determining the slope is still not unique. In Appendix A, some possible conditions are discussed. Based on the discussion given in Appendix A, we assume that the slope t of the curve at point 3 is determined by t = (Im4--m31m2+ ]m2-ml[m3)/(Im4-ma[+]m2-roll), (1)

where ml, m2, m3, and m4 are the slopes of line segments 12, 23, 34, and ~ , respectively. Under this condition, the slope t of the curve at point 3 depends only on the slopes of the four line segments and is independent of the interval widths. Under condition (1), t = m2 when ml = m 2 and m 3 ~ ~rt4, and t = m 3 when m~ = m4 and ml # m2, as desired. I t also follows from (1) that, when m2 = m3, t = mE = m3. Invariance of condition (1) under a linear scale transformation of the coordinate system is also obvious.

Journal of the Association for Computing Machinery, Vol. 17, No. 4, October 1970

A N e w Method of Interpolation and Smooth Curve F i t t i n g

591

When ml = m2 ~ ms = m4, the slope t is undefined under condition (1): t h e slope t can take any value between m2 and ms when ml approaches m2 and m4 approaches ms simultaneously. I t is a cornerstone of our new method t h a t t = mE when ml = m2 and, similarly, t = ms when m4 = m3, and these two rules conflict when ml = m 2 ~ ms = m4 ; therefore, no desired curve exists under condition (1) in this special case. (In order to give a definite unique result hi all cases, the slope t is equated to (m2 ~ ms) as a convention for this case in the computer programs, which are described in Section 4. This convention is also invariant under a linearscale transformation of the coordinate system.)
2.2. I N T E R P O L A T I O N BETWEEN A PAIR O F P O I N T S . W e t r y to express a portion of the curve between a pair of consecutive d a t a points in such a way t h a t the curve will pass through the two points and will have at the two points the slopes determined by the procedure described in Section 2.1. To do so, we shall use a polynomial because, as stated by Milne [6], "polynomials are simple in form, can be calculated by elementary operations, are free from singular points, are unrestricted as to range of values, m a y be differentiated or integrated without difficulty, and the coefficients to be determined enter linearly." Since we have four conditions for determining the polynomial for an interval between two points (Xl, Yl) and (x2, y2), i.e.

Y = Yl
y= y~

and
and

dy

~x = tl
dy ~x= th

at at

x-- xl,
x--Is,

where t~ and t2 are the slopes at the two points, a third-degree polynomial can be uniquely determined. Therefore, we assume t h a t the curve between a pair of points can be expressed by a polynomial of, at most, degree three. The polynomial, though uniquely determined, can be written in several ways. As an example we shall give the following form:
Y = Po + p l ( x -- xl) + p2(x -- x~)2 + p3(x -- x~) 3,

(2)

where
po = yl ,

(3)

pl = tl, P2 = [3(y2 -- y , ) / ( x 2 -- xl) -- 2tl -- t2]/(x2 -- xl), p3 -- It1 -~ t2 -- 2(y2 -- y l ) / ( x 2 -- xl)]/(x2 -- Xl)2.

(4) (5) (6)

2.3. ESTIMATION OF TWO MORE POINTS AT AN END POINT. At each end of the curve, two more points have to be estimated from the given points. We assume for this purpose t h a t the end point (x3, Ys) and two adjacent given points (x2, y2) and (xl, yl), together with two more points (x4, y4) and (x6, Ys) to be estimated, lie on a curve expressed by
y = go + gl(x Is) + g2(x -- x3)2,

(7)

where the g's are constants. Assuming t h a t

Journal of the Association for Computing Machinery, Vol. 17, No. 4, October 1970

592

HIROSHI AKIMA

we can determine the ordinates y4 and Y6, corresponding to x4 and xs, respectively, from (7). The results are
(Y5 --

y,)/(x~

- - x , ) - - (Y4 - - y ~ ) / ( x 4 - -

x3)
- - x3) - -

---- (Y4 - - y 3 ) / ( x , = (y3 - - y s ) / ( x 8

(Y3 -- y 2 ) / ( z a - - x2)
- - xx).

(9)

- - z2) - - (y2 - - y l ) / ( x 2

3.

Compar~on

With

Some Other Methods

Using a simple example taken from a study of waveform distortion in electronic circuits being conducted by the author, we compare our new method with four others. Assume t h a t the values of x a n d y at 11 points are as follows:
x 0

y =y(x)

1 2 3 4 5 6 7 8 9 10 10 10 10 10 10 10 10.5 15 50 60 85"

Knowing from the physical nature of the phenomena t h a t y ( x ) is a single-valued smooth function of x, we try to interpolate the values of y ( x ) and to fit a smooth curve to the given set of data points. First we apply the method of interpolation based on polynomials [4, 6]. This method is, perhaps, the one most often used. The result obtained by applying the tenth-degree polynomial is shown in Figure I(A). Second we use another well-known method based on the Fourier series (see [4, Sec. 9.3]). In applying this method to our data, we assume t h a t the whole range of x from 0 to 10 corresponds to one-half of the fundamental period from 0 to ~- and apply a series of cosine functions up to the tenth-order harmonic term. The result is shown in Figure 1 (B). The third method is based on a spline function [3]. The spline function of degree n is a piecewise function composed of a set of polynomials, each of degree n, at most, and applicable to successive intervals of the given data points. All the polynomials are determined as a set, so t h a t the function and its derivatives of order 1, 2, -.- , n - 1 are continuous in the whole range of x. Note that, although the spline function of degree three is somewhat similar to our method, no individual polynomial can be determined locally in the spline function. The result of applying the third-degree spline function is shown in Figure i(C). Next we try to apply the osculatory interpolation method [1, 5]. It is, like our method, based on a pieeewise function composed of a set of third-degree polynomials, each applicable to successive intervals of the given points, with the slopes at the given points locally determined. The only difference between this method and ours is the manner of determining the slopes at the given points. In the osculatory interpolation, the determination of the slope involves only three points, i.e. the data point in question as a center point and two neighboring data points. It is assumed t h a t the slope of the desired curve at the data point is equal to the slope at the same point of the curve of the second-degree polynomial passing through the three points involved. From the three data points given, (xl, yl), (x2, y2), and (x3, y~), the slope t at the center point (x2, Y2) is determined by t -- [(x2 -- Xl)2(y3 - - y2) + (x3 - - x , ) * ( y 2 - - y l ) ] /
[(x~ x,)~(x3 x2) + (x~ x~)~(x2 - - x2) +

x~)]

(10)

-- [(x3 -- x 2 ) m l + (x2 - - x l ) m , ] / [ ( x 3

(z2 -- xl)],

Journa| of the Association for Computing Machinery, Vol. 17, No. 4, October 1970

A New Method of Interpolation and Smooth Curve F~ing

593

/
(C) SPLINE {DEGREE 3} (D) OSCULATORY

__f
(E) NEW METHOD
(F} MANUAL

y.

{6 PEOPLE}

FiG. 1. Comparisonof several methods of smooth curve fitting. (Encircled points are given data points.) where ml and m2 are the slopes of line segments 1-'2 and 2-'3, respectively. The result of applying this method is shown in Figure I(D). Finally, we apply our new method, with the result shown in Figure i(E). In addition to these results of mathematical methods, the curve obtained manually is shown in Figure i(F); it is the average of curves drawn manually by six scientists and engineers. Compar~on of the curves in Figure i (A)-(F) indicates that the first two methods are definitely unsuitable for the example given. Although the curves obtained by the spline function method and the osculatory interpolation method, shown in Figure 1 (C) and (D), respectively, resemble the one obtained manually, shown in (F), they have maxima and minima that are absent in the manually drawn curve. The curve obtained by our new method, shown in Figure I(E), is closer than the other curves to the manually drawn curve in (F). In Figure 2 (A)-(F), we further compare our new method with the spline function method and the osculatory interpolation method with the same y values as
Journal of the Association for Computing Machinery, Vol. 17, No. 4, October 1970

594
(DEGREE 3)

HIROSHI AKIMA

(A~ SPLINE

(DEGREE 3)

(B) SPLINE

j
(C) OSCULATORY [O) OSCbtATORY (EJ NEW METHO0 (F) NEW METHOD

Fie. 2. Further comparison of three methods when the data points are given at unequal intervals. (Encircled points are given data points.) in Figure 1-(A)-(F), but with different x values. In Figure 2 (A), (C), and (E), the x values are 0, 1, 3, 4, 6, -.- , 13, and 15. T h e y are 0, 2, 3, 5, 6, . . . , 14, and 15 in Figure 2 (B), (D), and (F). T h e comparison indicates that, even when the data points are given at unequal intervals, our method performs as well as or better than the other two methods.

4. ComputerApplications
Our new method is further compared with the spline function method and the osculatory interpolation method from the standpoint of its computer applications. Two types of computer subroutines were programmed for each method, both in the CDC-3800 FORTRAN: one for interpolation and the other for smooth curve fitting. The interpolation subroutine interpolates, from a given set of coordinates of data points and for a given set of abscissa values of desired points, ordinate values of the

Journal of the Association for Computing Machinery, Vol. 17, No. 4, October 1970

A New Method of Interpolation and Smooth Curve Fitting

595

desired points. The curve fitting subroutine equally divides each interval between a pair of given points, interpolates an ordinate value for each dividing point, and generates a new set of points consisting of the given data points and the interpolated points. For the spline function method, the FORTRAN program given in [3] was used after slight modifications--modifications made to save both program length and computation time for fair comparison. 4.1. PROGRAMLENGTHS. The interpolation and smooth curve fitting subroutines based on our new method occupy 339 and 324 locations, respectively, compared with 245 and 276 locations occupied by the same subroutines based on the osculatory interpolation method. Program lengths of the same subroutines based on the spline function method depend on the maximum number of given data points Lm~x that can be processed by the subroutines. The subroutines, based on a thirddegree spline function, occupy 270 + 3Lm~x and 293 + 3Lmax locations. In each application, program length for our method is longer than that for the osculatory interpolation method. Comparison between our method and the spline function method depends on Lmax ; program lengths for these two methods are nearly equal for interpolation and smooth curve fitting at Lm~x = 25 and 10, respectively. 4.2. COMPUTATIONTIMES. Each subroutine was run many times on the CDC3800 computer, and the running time was measured by the internal clock. The subroutines based on the spline function method were run with the error tolerance in iterative solution of the equations for the second derivative of the spline function of 104 and the input data scaled in several ways, and the averages of the running times were taken. The results are shown in Tables I and II for interpolation and smooth curve fitting, respectively. For simplicity, we denote the number of given data points for both interpolation and curve fitting by L, the number of desired points for interpolation by N, and the number of divisions in each interval for smooth curve fitting by M. In each application, the time required by the osculatory interpolation is the shortest for a given combination of L and N or L and M. For interpolation, comparison between our new method and the spline function depends on the combination of L and N and also on the way the abscissa values of the desired points are given. Table I indicates that our new method almost always requires less time than the spline function method when the abscissa values of the desired points are given in an ascending order. For smooth curve fitting, the new method always requires less time than the spline function method.

5. Concluding Remarks We have described a new mathematical method of interpolation and smooth curve fitting. For proper application of our new method, the following remarks seem pertinent. (1) Since the curve obtained by our method passes through all the given points, the method is applicable only when the precise values of the coordinates of the data points are given. All experimental data have some errors in them, and unless the errors are negligible it is more appropriate to smooth the data, i.e. to fit a curve approximating the data appropriately, than to fit a curve passing through all the points.
Journal of the Association for Computing Machinery, Vol. 17, No. 4, October 1970

596

HIROSHI AKIMA

TABLE I. COMPARISONOF COMPUTATION TIMES OF INTERPOLATION BY VARIOUS METHODS. (L = NUMBER OF GIVEN DATA POINTS AND N = NUMBER OF DESIRED POINTS.)
L N
Computation times (msec)

Spllne

Osculatory

New method

10

1 10" 10"* 100" 100"* 1 10" 10"* 100" 100"*

20

40

10" 10"* 100" 100"* 1 10" 10"* 100" 100"*

100

10" 10"* 100" 100"*

3.0 4.6 4.7 17 21 7.2 9.2 9.6 24 30 20 23 23 42 48 46 49 50 75 81 120 125 125 160 160

0.8 2.0 2.3 12 15 0.8 3.0 3.6 17 27 0.9 3.5 4.3 22 33 1.0 4.6 5.4 32 47 1.0 5.3 5.8 41 53

1.0 2.6 3.3 16 26 1.1 3.9 5.2 21 46 1.2 4.4 5.8 27 55 1.2 6.3 7.6 36 67 1.2 7.0 8.0 49 74

* When abscissa values of the desired points are given in ascending order. ** When those values are given in random order. (2) A s is t r u e for a n y m e t h o d of i n t e r p o l a t i o n , t h e a c c u r a c y of t h e interpolation c a n n o t be g u a r a n t e e d , unless t h e m e t h o d in q u e s t i o n h a s been c h e c k e d in advance a g a i n s t precise v a l u e s o r a f u n c t i o n a l form. (3) O u r m e t h o d y i e l d s a s m o o t h , n a t u r a l - l o o k i n g c u r v e a n d is t h e r e f o r e useful in cases w h e r e m a n u a l , b u t tedious, c u r v e f i t t i n g will do in principle. (4) T h e r e s u l t a n t c u r v e of o u r m e t h o d is i n v a r i a n t u n d e r a linear-scale transf o r m a t i o n of t h e c o o r d i n a t e s y s t e m . I n o t h e r words, different scalings of the coo r d i n a t e s r e s u l t in e q u i v a l e n t curves. (5) O u r m e t h o d is nonlinear. I n o t h e r words, if y~ = y / ~ y t~ for all i, t h e interpol a t e d v a l u e s d o not, in general, s a t i s f y y ( x ) = y ' ( x ) ~ y'~(x). (6) O u r m e t h o d gives e x a c t r e s u l t s w h e n y is a s e c o n d - d e g r e e p o l y n o m i a l of x, p r o v i d e d t h e abscissas of t h e d a t a p o i n t s a r e e q u a l l y spaced. (7) O u r m e t h o d r e q u i r e s o n l y s t r a i g h t f o r w a r d procedures, n o t i t e r a t i v e solutions of e q u a t i o n s w i t h p r e a s s i g n e d e r r o r tolerances, w h i c h are r e q u i r e d b y some methods. N o p r o b l e m concerning c o m p u t a t i o n a l s t a b i l i t y o r convergence exists in applicat i o n of o u r m e t h o d . (8) O u r m e t h o d can b e i m p l e m e n t e d as c o m p u t e r s u b r o u t i n e s w i t h reasonable

Journal of the Associationfor ComputingMachinery,Vol. 17, No. 4, October1970

A New Method of Interpolation and Smooth Curve Fitting

597

COMPARISON OF COMPUTATION TIMES OF SMOOTH

FITTING

(L

NUMBER

OF

GIVEN DATA POINTS AND M NUMBER OF DIVISIONS IN EACH INTERVAL BETWEEN A PAIR OF SUCCESSIVE DATA POINTS.)
L M

Computatlon TimeJ (reset)

10

20

40

100

2 10 100 2 10 100 2 10 100 2 10 100 2 10 100

3.3 4.9 21 8.0 12 50 22 29 110 50 64 230 130 170 590

1.5 3.0 20 2.8 6.1 44 5.2 12 93 10 25 190 25 63 480

2.0 3.6 21 3.7 6.9 46 7.0 14 95 13 28 200 33 71 500

3 C A

O /9

FIG. 3. A condition for determining the slope of the t a n g e n t to the curve, previously proposed [2]

program length. There is no drawback in computation time compared with typical existing mathematical methods. (9) Our method, as reported on in the text of this paper, is applicable only to a single-valued function, but the basic idea of our new method is also applicable to the case of a multiple-valued function. A modified method for this case is outlined in Appendix B.

Appendix A. Curve

A Supplementary Note on the Determination of the Slope of the

With five data points 1, 2, 3, 4, and 5, as shown in Figure 3, the author [2] previJournal of the Association for Computing Machinery, Vol. 17, No. 4, Ootober 1970

598

H I R O S H I AKIMA

ously proposed a condition t h a t the slope of line segment ~ , at point 3, is determined by

tangent to the curve (11)

I~/~Xl

= 14D/DB I,

where the point of intersection of the two straight lines extended from line segments 1-2 and 3--4is denoted by A, the same point corresponding to line segments ~ and 45 by B, and the points of intersection of the tangent with two straight lines extended from 12 and 45 are denoted by C and D, respectively. The ratio of two line segments, which is used when the two line segments are on a straight line, is used here to symbolically represent the ratio of the lengths of the two line segments with a plus or minus siga depending on whether these two line segments have the same sense or not. An analytical expression of the slope of CD under condition (11) will be derived in the next paragraph. Let the coordinates of points 1, 2, 3, 4, 5, A, B, C, and D in Figure 3 be denoted by (x, y) with subscripts 1, 2, 3, 4, 5, a, b, c, and d, respectively. For simplicity, we define ai = xi+l -- xi
b,

(i = 1, 2, 3, 4), (i = 1, 2, 3, 4).

(12) (13)

= y~+~ - y~

We denote the slopes of 12, 23, 34, 45, and CD by m l , m2, m3, m4, and t, respectively, i.e.
m~ = (y~+l y~)/(x~+~ t = x~) = b~/a~ -Xc).

(i = 1, 2, 3, 4),

(14) (15)

(Yd --

y~)/(Xd

Then, it is clear from Figure 3 t h a t the following equations should hold:

(16) (17) (18) (19)

-- Xb) ---X~)

(Y4 - b3/a3, b2/a2, (Yd --

(ya -(Yb --

X3)
Xc)

(Y3 --

y3)/(Xd

--

X3) = t.

(20)

Using (12) and (13), we have, from (16), [(y8

(21)

Eliminating (Y3 -- Ya) from (18) and (21), we obtain

(22)

Following the same procedure as in obtaining (22) from (16) and (18), we obtain
(xb xa)/a2 x3 Xd -xc = = (a3b4 - (alb2 -a,b3)/(a~b4 a~bl)/(a~t a463)/(b4 --a,b~), bl), a4t).

(23) (24) (25)

X3 =

(a364 - -

Journal of the Association for Computing Machinery, Vol. 17, No. 4, October 1970

A N e w M e t h o d of I n t e r p o l a t i o n a n d S m o o t h Curve F i t t i n g

599

Since condition (11) can be w r i t t e n as [[(x3 -- x~) -- a2]/[(xa -- x~) -- (xa -- x~)]] = I[aa-- (xd-za)]/[(Xd--

Xa) -- ( X b - - X3)][,

(26)

we obtain, by eliminating (xa -- Xa), (Xb -- Xa), (Xa -- X), and (xe -- xa) from (22)(26), a quadratic equation of the form [ S~2S2t ] (eat -- ba) 2 = [S~aSa, [ (a2t where
S~i = a ~ b j asb~ (i ~ j ) .

52) ~,

(27)

(28)

From the r e q u i r e m e n t t h a t points 2 and 4 m u s t lie on t h e same side of the t a n g e n t to the curve at point 3, it follows t h a t
(a2t -- b2)(aat -- ba) < 0.

(29)

Therefore, we h a v e I S12Z24 ](ba -- ant) = I SlaZa, Ii(a2t Solving this equation, we obtain t ---- (w2b2 + wab3)/(w2a~ + waaa), where
w~ = [ S~Sa4 [, wa = ] S12S~4 1.

52).

(30)

(31) (32) (33)

This is an analytical expression of the slope of C D determined under condition (11). The relations (31)-(33) can also be expressed as t = (w2'm2 -k w 3 ' m a ) / ( w 2 ' + wa'), where w2' = (sign of as) I (m3 -- m l ) ( m 4 -- ma) 1, wa' = (sign of aa) I (m2 -- ml)(m4 -- m2) [~. (35) (36) (34)

Note that, in the case of a single-valued function, we can always m a k e a2 and aa positive and simplify (35) and (36). Since condition (11) is a purely geometrical one, it is i n v a r i a n t under a linear transformation of the coordinate system, which includes a linear scale transformation and a rotation. I t follows from (34)-(36) and also from the geometrical construction dictated b y (11), t h a t the slope t depends only on the slopes of four secants, i.e. on the quantities m l , m 2 , m 3 , and m , , and is independent of the interval widths. I t is also clear from (34)-(36) t h a t ml = m2, m3 ~ m l , and m4 ~ m3 implies t = ml = m2, and similarly t h a t m3 = m4, ml ~ m2, and m4 ~ ms implies t = ma = m4 : these are highly desirable properties. W h e n m l = m 2 ---- m a or when m2 = m3 = m4, t is undefined b y (34); however, this difficulty can easily be settled by taldng t = m2 = ma for such cases as a n a t u r a l and logical extension. However, despite its desirable properties, condition (11) has a serious drawback.

Journal of the Association for Computing Machinery, Vol. 17, No. 4, October 1970

600

HIROSHI

AKIMA

I t follows from (34)-(36) t h a t t = m2 when m2 = m4, m3 ~ ml, and m4 ~ m3, and similarly t h a t t = m3 when m3 = ml, m2 ~ ml, and m~ ~ m2. These properties are by no means desirable or to be expected. The simplest way to eliminate these undesirable properties for the case of a single-valued function is, perhaps, to modify weighting coefficients w2' and w3' in (34) to read we' = Ira4 , m3 I,
m ' = l m2 mll,

(37)
(38)

instead of (35) and (36), respectively. By this modification, most of the desirable properties of condition (11) are retained. The only exception is t h a t the property of invariance under a rotation of the coordinate system is lost. But the requirement for invariance under a rotation is immaterial in the case of interpolation of a singlevalued function. Equation (34) with w2' and wa' defined by (37) and (38), respectively, is used in the text of this paper as the equation for determining the slope. This is equivalently given in the text as condition (1). In the case of a multiple-valued function, on the other hand, the property of invariance under a rotation of the coordinate system is essential. While retaining this property, the undesirable properties of (11) can be eliminated by modifying w2 and w8 in (31) to read w2 = ]$34 ],
w, = I s , 2 I,

(39)
(40)

instead of (32) and (33), respectively. Invarianee under a rotation by using (31) with these modified w2 and w3 can be realized from the fact t h a t [ S~ [ is an invariant quantity representing twice the area of the triangle bounded by the vectors (a, b~) and (a, bi), and t h a t (31) can be written in the form

&,/,s,,

IS=l/I,S~ I,
atb2,

(41)

if at and bt are chosen so t h a t t = bt/at and Sst and St3 are defined by S2t = a 2 b t -

(42) (43)

St3 = atba - a3bt,

respectively, as S~ is defined by (28). These modified w coefficients in (39) and (40) correspond to
W2' = a2 l a3a4(ma -- m3) l, w d = a8 [ alaffm2 ml)1.

(44) (45)

I t follows from (44) and (45) t h a t another property of dependence of t on m l , m2, m3, and m4 only is lost, but this property is not important in the case of a multiplevalued function. Equation (31) with the weighting coefficients defined by (39) and (40) is used in Appendix B.
Journal of the Association for Computing Machinery, ol. 17, No. 4, October 1970

A N e w Method of Interpolation and Smooth Curve Fitting Appendix B. Function

601

A n Outline of the Method of Interpolation for a Multiple-Valued

As described in Appendix A, (31) with (39) and (40) can be used in principle to determine the slope of the curve in the case of a multiple-valued function. To present the method in a form t h a t is ready to be implemented in a computer program, however, it is necessary to rewrite all the equations without a ratio of the increment in y to t h a t in x, such as m~ or t in Appendix A. Since (39) and (40) are given in terms of the increments in the x and y directions, we have only to rewrite and eliminate t from (31). This can be done b y the use of cos 0 and sin 0 instead of t = tan 0, where 0 is the angle of the tangent to the curve measured from the x axis. The results are cos 0 = a0(a02 -~- bd) ~ , sin 0 = bo(ao2 + bo2)-. where
ao
= w2a2

(46) (47)

"~- w3a3,

(48) (49)

bo = ff)2b2 ~- robs,
w2 = [ $34 I = [ a3b, -- a453

I,

(50)
(51) (52) (53)

w3 = [ $12 I = [ alb2 -- asbl I, a, = x~+l -- x~ (i = 1, 2, 3, 4),

b~ = y~+l - Y~

(i = 1, 2, 3, 4).

We assume t h a t the curve between a pair of d a t a points (xl, yl) and (x2, y~) can be expressed b y
x = po + plz ~- p2z 2 ~ p3z3, Y = qo -~ qlz -~ q2z~ -~ q3zS, (54)

(55)

where the p ' s and q's are constants, and z is a p a r a m e t e r t h a t varies from 0 to 1 as the curve is traversed from (xl, yl) to (x2, y2). Since the coordinates of the two points (Xl, yl) and (x2, y~), as well as the direction of the curve (cos 0~, sin 01) and (cos 02, sin 02) at these points, are given, we further assume t h a t x and y satisfy the condition
X =Xl, X ~X2,

y = yl, Y = y2,

dx/dz = rcosOj, dx/dz = rcosO~,

and and

dy/dz = rsin01 dy/dz = rsin02

at at

z = 0, z = 1,

where r = [(x2 -- xl) 2 + (Y2 -- Yl)2]~. (56)

From these conditions we can uniquely determine the p and q constants. The results are p0 = x l , p, = r cos 01, (57) (58)

Journal of the Association for Computing Machinery, Vol. 17, No. 4, October 1970

602 p2 = 3 (x2 -- xl) -- r (cos 02 + 2 cos 01),

HIROSHI AKIMA

(59)

(60)
(61)

qo = y l , ql = r s i n 0 1 , q2 = 3 (Y2 -- Yl) -- r (sin 02 --]- 2 sin t~l), q3 = - - 2 (y~ -- y~) + r (sin ~2 -~- sin t~).

(62) (63) (64)

E x c e p t for the case of a closed curve, e s t i m a t i o n of two more points from the data points is required at each end of the curve. W e assume for this purpose t h a t the end point (x3, Y3) and t w o a d j a c e n t given points (x~, y~) and (xl, yl), t o g e t h e r with two more points (x4, y4) and (xs, ys) to be estimated, lie on a curve expressed b y
x = go glz ~ g2z 2,

(65) (66)

y = ho ..~ hlz ~ h2z 2,

where t h e g's and h's are constants and z is a p a r a m e t e r . Assuming t h a t x = x~ and Y = Yl at z = i (i = 1 , 2 , 3 , 4 , 5 ) ,

we can determine the g and h constants and, consequently, also the coordinates (x4, y4) and (x~, Ys). T h e results can be expressed as ( x s - - x4) -- ( x 4 - - x3) = ( x 4 - - x3) -- ( x 3 - - x~) = ( x 3 - - x2) -- ( x 2 - - xl), ( Y s - - y4) -- ( y , - - Y3) = ( y t - - Y3) -- (Y3-- Y~) = (Y3-- Y2) -- ( y 2 - - yl).
(67)

(68)

T h e interpolated curve is i n v a r i a n t u n d e r a r o t a t i o n of the coordinate system but v a r i a n t under a linear-scale t r a n s f o r m a t i o n . Therefore, b o t h t h e abscissa and the ordinate should be scaled with their respective units h~ving an equal actual length on t h e graph. ACKNOWLEDGMENTS. T h e a u t h o r expresses his deep appreciation to E. L. Crow, R. K. Rosich, J. S. W a s h b u r n , and R. E. Wilkerson of the I n s t i t u t e for Telecomm u n i c a t i o n Sciences, E S S A R e s e a r c h Laboratories, for their helpful discussions.
REFERENCES ACKLAND, T . G . On osculatory interpolation, where the given values of the function are at unequal intervals. J. Inst. Actuar. ~9 (1915), 369-375. 2. AKIMA, H. A method of smooth curve fitting. ESSA Tech. Rep. E R L 101-ITS 73. U S Government Printing Office, Washington, D. C., Jan. 1969. 3. GREVILLE, T. N . E . Spline functions, interpolation, and numerical quadrature. In Mathematical Methods for Digital Computers, Vol. 2, A. Ralston and H. S. Wilf (Eds.), Wiley, New York, 1967, Ch. 8. Introduction to Numerical Analysis. McGraw-Hill, New York, 1956, 4, HILDEBRAND, F . B . Ch. 2, 3, 4, and 9. 5. KARUP, J. On a new mechanical method of graduation. In Transactions of the Second International Actuarial Congress. C. and E. Layton, London, 1899, pp. 78-109. 6. MILNE, W. E. Numerical Calculus. Princeton U. Press, Princeton, N. J., 1949, Ch. III.
1.

RECEIVED DECEMBER, 1969; REVISED APRIL, 1970

Journal of the Association for Computing Maohlnery, Vol. 17, No. 4, October 1970