- Email: [email protected]
Newton-Everett interpolation of continuous functions
Computer Physics Communications 27 (1982) 73—78 North-Holland Publishing Company
73
NEWTON-EVERETF INTERPOLATION OF CONTINUOUS FUNCTIONS J.A. HERNANDO Depto. de Fisica, Comisión Nacional de Energia Ató,nica, Av. del Libertador 8250, 1429 Buenos Aires, Argentina Received 3 November 1981
PROGRAM SUMMARY Title of program: FINT
interpolation routine that employs the finite difference for-
Catalogue number: AARS
mulae of Newton and Everett [I]. The routine is designed to work with a continuous function defined on an evenly spaced set of points. The main purpose is to do the calculations in a quick and accurate form. The
Program obtainable from: CPC Program Library, Queen’s University of Belfast, N. Ireland (see application form in this issue) Computer: IBM 370/158; Installation: Comisiôn Nacional de Energia Atómica, Av. del Libertador 8250, 1429 Buenos Aires, Argentina Operating system: OS/VSI or VM370/CMS Programming language used: FORTRAN IV
calculation is divided into two parts: i) a calculation which depends exclusively on the points where the function is given and on the accuracy requested; ii) a calculation dependent on the point where one wishes to interpolate. The first calculations are done only once and are then
High speed storage required: 234 Kbytes
available for the following calculations.
No. of bits in a word: 32
Restrictions on the complexity of the program In the program it is assumed that the set of points where the function is defined is less than or equal to 200. This can be
Overlay structure: none No. of magnetic tapes required: none No. of cards in combined program and test deck: 673
easily changed. If the function to be interpolated has a pathological be-
Card punching code: EBCDIC
haviour, the results are strongly dependent on the size of the separation between functional values, so special care must be
Keywords: general purpose, Newton interpolation, Everett interpolation, finite differences, continuous approximation
taken with this point (see section 3). Typical running time
Nature ofphysical problem Interpolation needs are deeply rooted in almost all branches of physics because of the uncommonness of analytical solutions and the discrete nature of physical measurements.
The average running time depends on the number of interpolation points calculated. However, it can be estimated to be of the order of 3 times that needed for a double-precision exponential (see section 3).
Method of solution If the calculations one is doing are very time consuming, it would be convenient to minimize these calculations and to interpolate whenever possible. In this paper we describe an
Reference [I] Z. Kopal, Numerical analysis (Chapman and Hall, London,
001 0-4655/82/0000—0000/$02.75 © 1982 North-Holland
74
iA. Hernando
/
Newton —Everett interpolation of continuous functions
LONG WRITE-UP
1. Introduction
2. The interpolation formalism [1,3]
In the days of hand calculation the need of to get easy access to numerical values of the most diverse functions encouraged the construction of tables as accurate as possible within a reasonable size. Therefore, interpolation formalism was a rapidly developing branch of numerical analysis. However, the explosive development of computer technology led to the rather general belief that the days of interpolation calculus were past (see ref. [1], page 40). This is true for the large class of problems that involve the evaluation of some “good” function (good in the sense of having at our disposal a fast algorithm to evaluate it). But there are many problems that do not fit in this scheme; an important example is the solution of a non-linear integral equation such as those appearing in the statistical theory of liquids [2]. In such a case the first decision that must be taken is which integration method to employ. Roughly speaking, we have two kinds of methods: i) those in which the function to be integrated is calculated on a set of evenly spaced points (Newton—Cotes and Simp-
The function to be considered is defined in a closed interval. We know the values on a set of evenly spaced points and we need to interpolate the function at a generic point of the interval. This point can be written as ~ + h ~ ~ 1 (1) p where x~is the s th point at which we know the function, h is the separation between the functional points and p indicates the fractional separation from x~.The Newton-forward interpolation formula to order n in the differences is:
son s methods are the most widely known); ii) those in which the function to be integrated is calculated on a set of irrational abscissas (e.g. Gauss and Tchebysheff’s methods). Integration methods of the second group are, by far, faster and more accurate than methods of the first group. Therefore, one clearly sees the convenience of implementing them in any iterative search for the solution to the integral equation. As a consequence, it is crucial to have an interpolation routine for evaluating the function to be integrated at the irrational abscissas needed. This routine must be such as not to spoil the advantages of these integration methods. In this article we describe an interpolation routine that employs the Newtonforward, Newton-backward and Everett interpolation formulae. These formulae are displayed in section 2. In section 3 we describe the routine. In section 4 some general remarks are included.
.~
‘
f(x, +ph)
~
+
(P) k
~
k=
O~p~ 1,
~kj
(2)
the Newton-backward formula is
—
I p+k
~
fc~x +ph)—f+
~ i k= I
—
1
/i vkf
k
~‘
I ~p ~ 0
(3)
and the Everett interpolation formula is (n even) f(x
h)
+ s
+
P
k—I
f/I —p + k \ 2k + I
—
+
2k
)~
(~2kp ++k1
2kf
6
J
0~
~ I P
‘
(4) In these eqs. J~=f(x~); the combinatorial factors are defined by (P) = p (p I)... (p k + 1)/k! k and the finite differences are —
(5)
—
I.
forward
~f
=
~ (~ 1)~(k\f~~
~J/
j0 k
backward
vkf
central
(6a)
~‘
(~)
~ (— 1)~ f— 1=o 12k
,
(6b)
fs+k—r
(6c)
2k 82kf
~ (—l)~~ ~ j=0
J.A. Hernando
/ Newton
—
Everett interpolation of continuous functions
It must be noticed that the forward difference of order k contains all the function values between f, and fs+k’ the backward difference of order k contains the values between f~k and f~,and the central difference of order 2k, centred in J, contains the values between £—k and J.+k. Therefore, it is clear that the Everett interpolation is applicable only in the middle part of the interval and the Newton formulae are specially suited for the ends of the interval, The performance of these interpolation formulae is strongly dependent on the function to be interpolated, on the size of h and on the accuracy needed. This will be discussed in the next section. However, an important point must be emphasized: the tabular differences are a sensitive detector of roundoff and/or tabular errors [1,3]. In this discussion we consider only roundoff errors. When the regular flow of functional values is perturbed by random errors in some decimal place, it manifests itself in the appearance of an oscillating pattern in some tabular difference. The order of this finite difference tells us the accuracy of the functional values and the greatest finite difference that it makes sense to take into account.
3. The program The program is a double-precision function named FINT(XINT). The routine employs the Everett interpolation formula whenever possible and Newton formulae at the interval ends. If we look at the formulae (2), (3) and (4) we see that (for fixed n) we can rearrange them in a power series in p. In such a case each summand is a product of two factors: one solely dependent on the differences (i.e. on the subdivision of the interval and on the order of the interpolation) and the other being some power of p. Therefore, the factor dependent on the differences is evaluated only once in the first call to FINT. In all the subsequent calls these quantities are available. The distinction between these two cases is made with the logical command DIFF. When DIFF is true (first call), the quantities dependent on the differences are evaluated and control is returned to the main program without calculating any interpolated
75
value. In this case XINT is a dummy variable. When DIFF is false, the differences are supposed known and the program jumps to the direct evaluation of p and of the interpolated function at this point. It must be noticed that this methodology is efficient if we need to evaluate several interpolated values (whose number is at least equal to the number of known functional values NDIM, see below). In any other case it is more time consuming, than the more conventional method that evaluates exclusively the differences needed for each point and sums the series as it stands in (2), (3) or (4). There are other five logical commands: TWO, FOUR, SIX, EIGHT, TEN which direct the program to evaluate the function to the requested finite difference order. The last logical variable, WARN, is true when an error condition has been detected; in such a case the user must decide which action should be taken by the main program. All these logical commands are in the labelled common INT3. The functional values and the abscissas to which they correspond are in the vectors XCALC, FCALC. These two vectors and the separation between abscissas (H) are in the labelled common INTl. The dimension of these vectors (NDIM) and the integer LP are in the labelled common INT2. LP is an integer supplied by the user to indicate the output device where error messages should be sent. The error messages are of three kinds. The first one appears in the case that the five logical cornmands TWO,.. .,TEN are all false. The second message is for the case that the number of functional values (NDIM) does not suffice to evaluate the interpolation to the requested finite differences (NDIM must be greater than 2, 4, 6, 8, 10 for second, fourth, sixth, eighth and tenth interpolation order, respectively). The third error message is for the case that the point XINT, at which we want to evaluate the interpolated function, lies outside the interval covered by XCALC. There are several points worthy of mention with respect to error conditions. A first point concerns the change of interpolation orders (e.g. from TWO to FOUR). In this case all the quantities dependent on the differences should be recalculated;
76
J.A. Hernando
/ Newton —Everett interpolation of continuous functions
Table I CPU times and relative errors for the various options CPU time needed when: (for the Units, see text)
Relative error in the interpolated function when: _________________________________________________________________
f=exp(—x) cos(2~rx) DIFF=.T.
(T
1)
TWO FOUR SIX EIGHT TEN
II 24 48 78 120
___________________
___________________
(T,)
H=0.03
11=0.08
11=0.08
11=0.3
38 ‘U 43 46 48
tO ~ I0~ tO 8
10 iO~
to
to
lO~
~
IOb0
tO~
tO ~
tO
2
10 6 0 ~ 10 ~
10
~o
to
‘5
therefore, DIFF must be set true. A second point is that the routine does not check whether the vector XCALC is an ordered and equispaced vector; it assumes that to be true. If this is not true, the interpolation formulae (2), (3) and (4) do not apply. Another point concerns the possibility of having two commands true; e.g. FOUR and EIGHT. The routine does not detect this condition and performs the calculations ordered by the greatest difference command (EIGHT in this case). This superposition can be a source of errors, only some of which are covered by the checks included (e.g. that NDIM does not suffice for eighth order interpolation). With respect to the LP parameter: there is no default value, it must always be supplied by the user, The program is rather long and expensive in storage because our primary aim was to obtain an accurate and fast routine at the cost of increasing the storage requirements. Table 1 gives the CPU times for the various options that the program has built in. They are measured in units of the time needed to do a double precision multiplication, so that it can be estimated for another computer. If T1 is the time spent when DIFF is true (first call) and T2 when DIFF is false; the average time spent by interpolation in N calls to FINT is T
f=exp(—x)
DIFF=.F.
(NDIM*T 1 +N*T2 )/N.
In our system the time needed for a double precision exponential is 25 in these units, so it is seen that this routine is rather fast. In table 1 we also show the relative accuracy of
‘
to
the interpolation for two types of functions with two values of H. It is seen that we need relatively few values of FCALC to obtain a satisfactory accuracy.
4. Conclusions Very recently [4] a program has been published that evaluates the Tchebysheff series approximation to continuous functions. Tchebysheff series have a very attractive feature: they minimize the maximum error of the approximation (minimax approximation). Also interesting is the fact that to evaluate the Tchebysheff series one does not need an evenly spaced set of functional values. On the other hand, their coefficients must be determined by numerical integration and this implies both computation time and the dilemma of which integration method must be used. With the trapezoidal methods we lose accuracy and time and with the other methods an interpolation formalism in the spirit of this paper is needed. Therefore, it seems to us that these two routines can complement one another. To summarize, the main purpose of this paper was to present a fast and accurate interpolation routine. With an behaviour adequate choice of H (strongly dependent on the of the function, as can be seen in table I), it is possible to interpolate at practically the full accuracy of the given function. The economy in computation time increases with the number of interpolations calculated.
J.A. Hernando / Newton — Everett interpolation of continuous functions
References [I] A. Ralston, A first course in numerical analysis (McGrawHill, New York, 1965). [21 J.P. Hansen and I.R. McDonald, Theory of simple liquids (Academic Press, New York, 1976).
77
[3] Z. Kopal, Numerical analysis (Chapman & Hall, London, 1961). [41M.A. Christie and K.J.M. Moriarty, Comput. Phys. Commun. 23 (1981) 145.
78
J.A. Hernando
/
Newton — Everett interpolation of continuous functions
TEST RUN OUTPUT NEWTON—EVERETT INTERPOLATION FORMULAE FROM SECOND TO TENTII TABULAR DIFFLRENCF.S ORDER FUNCTION TO BE INTERPOLATED, FEXP(—X)COS(2PIX) 01FF ORDER TWO TWO TWO TWO
DIFFERENCE BETWEEN X VALUES, 9= 3.000—02
TWO
X VALUE 2.30000—02 l.4600D—0l 2.69i0D—01 3.92000—01 5.l5OiD—Ol 6.3800D—01 7.61000—01 8.8400D—01
FUNCTION VALUE= 9.670756160—01 5.25347650D—01 —9.100729170—02 —5.260101800—il —5.948488660—01 —3.418707650—01 3.226477460—02 3.001684980—il
INTERPOLATION= 9.669253700—01 5.25346730n—0l —9.100553800—02 —5.26005449D—01 —5.948419290—il —3. :18702900—il 3.2~560014D—02 3.081127220—01
RELATIVE ERROR= —l.554D—04 —1.7510—06 —1.9260—05 —8.9930—06 —1.1660—05 —1.3800—06 —2.719D—04 —1.810D—04
FOUR FOUR FOUR FOUR FOUR FOUR FOUR FOUR
2.30000—02 1.46000—01 2.6900D—01 3.92i0D—il 5.1500D—il 6.38000—01 7.61000—01 0.8400D—01
9.670756160—il 5.253476500—01 —9.10072917D—02 —5.260101800—01 —5.940488660—01 —3.41070765D—01 3.226477460—02 3.081684980—01
9.670784790—01 5.253476700—il —9.10072746D—02 —5.260101470—01 —5.948488330—01 —3.418707900—01 3.226468810—02 3.08168978D—01
2.96112—06 3.0210—08 —1.073D—i7 —6.229D—08 —5.5740—08 7.5210—08 —2.6790—06 1.5570—06
SIX SIX SIX SIX SIX SIX SIX SIX
2.30000—02 1.46000—01 2.69000—01 3.92000—01 5.lSOiD—0l 6.38000—01 7.61000—01 8.84000—il
9.670756160—01 5.253476500—01 —9.100729170—02 —5.260101800—01 —5.948488660—01 —3.418707650—01 3.226477460—02 3.081684980—01
9.670755570—01 5.253476500—01 —9.100729150—02 —5.26010l8iD—0l —5.948488660—01 —3.418707650—01 3.226477380—02 3.081684960—01
—6.1060—08 6.7610—10 —1.6880—09 —4.0980—10 —1.5280—10 1.2110—09 —2.4120—08 —6.6200—09
FlOE EIGH EIOH EIGH EIGH EIGH EIGH 12109
2.30000—02 1.46000—01 2.69000—01 3.92000—01 5. l500D—01 6.38000—01 7.61000—01 8.84000—01
9.670756160—01 5.2534765iD—01 —9.100729170—02 —5.260101800—01 —5.948488660—01 —3.418707650—01 3.226477460—02 3.081684980—01
9.670756170—01 5.253476500—01 —9.100729170—02 —5.260101800—01 —5.948408660—01 —3.418707650—01 3.226477460—02 3.081684970—01
1.1590—09 8.015D—12 —1.4150—11 —2.3620—12 1.2010—12 1.391D—1l —2.0090—10 —4.9520—10
TEN TEN TEN TEN TEN TEN TEN TEN
2.30000—02 1.46000—01 2.69000—01 3.92000—01 5.15000—01 6.38000—01 7.61000—01 8.84000—01
9.670756160—0]. 5.253476500—il —9.100729170—02 —5.260101800—01 —5.948488660—01 —3.418707650—01 3.226477460—02 3.081684980—01
9.670756160—01 5.253476500—il —9.l0i72917D—02 —5.261101800—01 —5.948488660—01 —3.41870765D—il. 3.226477450—02 3.081684980—01
—1.6760—11 5.8910—12 —1.1120—13 —9.89412—15 2.935D—14 1.3780—13 —5.3890—li 3.095D—ll
TWO
TWO TWO
73
NEWTON-EVERETF INTERPOLATION OF CONTINUOUS FUNCTIONS J.A. HERNANDO Depto. de Fisica, Comisión Nacional de Energia Ató,nica, Av. del Libertador 8250, 1429 Buenos Aires, Argentina Received 3 November 1981
PROGRAM SUMMARY Title of program: FINT
interpolation routine that employs the finite difference for-
Catalogue number: AARS
mulae of Newton and Everett [I]. The routine is designed to work with a continuous function defined on an evenly spaced set of points. The main purpose is to do the calculations in a quick and accurate form. The
Program obtainable from: CPC Program Library, Queen’s University of Belfast, N. Ireland (see application form in this issue) Computer: IBM 370/158; Installation: Comisiôn Nacional de Energia Atómica, Av. del Libertador 8250, 1429 Buenos Aires, Argentina Operating system: OS/VSI or VM370/CMS Programming language used: FORTRAN IV
calculation is divided into two parts: i) a calculation which depends exclusively on the points where the function is given and on the accuracy requested; ii) a calculation dependent on the point where one wishes to interpolate. The first calculations are done only once and are then
High speed storage required: 234 Kbytes
available for the following calculations.
No. of bits in a word: 32
Restrictions on the complexity of the program In the program it is assumed that the set of points where the function is defined is less than or equal to 200. This can be
Overlay structure: none No. of magnetic tapes required: none No. of cards in combined program and test deck: 673
easily changed. If the function to be interpolated has a pathological be-
Card punching code: EBCDIC
haviour, the results are strongly dependent on the size of the separation between functional values, so special care must be
Keywords: general purpose, Newton interpolation, Everett interpolation, finite differences, continuous approximation
taken with this point (see section 3). Typical running time
Nature ofphysical problem Interpolation needs are deeply rooted in almost all branches of physics because of the uncommonness of analytical solutions and the discrete nature of physical measurements.
The average running time depends on the number of interpolation points calculated. However, it can be estimated to be of the order of 3 times that needed for a double-precision exponential (see section 3).
Method of solution If the calculations one is doing are very time consuming, it would be convenient to minimize these calculations and to interpolate whenever possible. In this paper we describe an
Reference [I] Z. Kopal, Numerical analysis (Chapman and Hall, London,
001 0-4655/82/0000—0000/$02.75 © 1982 North-Holland
74
iA. Hernando
/
Newton —Everett interpolation of continuous functions
LONG WRITE-UP
1. Introduction
2. The interpolation formalism [1,3]
In the days of hand calculation the need of to get easy access to numerical values of the most diverse functions encouraged the construction of tables as accurate as possible within a reasonable size. Therefore, interpolation formalism was a rapidly developing branch of numerical analysis. However, the explosive development of computer technology led to the rather general belief that the days of interpolation calculus were past (see ref. [1], page 40). This is true for the large class of problems that involve the evaluation of some “good” function (good in the sense of having at our disposal a fast algorithm to evaluate it). But there are many problems that do not fit in this scheme; an important example is the solution of a non-linear integral equation such as those appearing in the statistical theory of liquids [2]. In such a case the first decision that must be taken is which integration method to employ. Roughly speaking, we have two kinds of methods: i) those in which the function to be integrated is calculated on a set of evenly spaced points (Newton—Cotes and Simp-
The function to be considered is defined in a closed interval. We know the values on a set of evenly spaced points and we need to interpolate the function at a generic point of the interval. This point can be written as ~ + h ~ ~ 1 (1) p where x~is the s th point at which we know the function, h is the separation between the functional points and p indicates the fractional separation from x~.The Newton-forward interpolation formula to order n in the differences is:
son s methods are the most widely known); ii) those in which the function to be integrated is calculated on a set of irrational abscissas (e.g. Gauss and Tchebysheff’s methods). Integration methods of the second group are, by far, faster and more accurate than methods of the first group. Therefore, one clearly sees the convenience of implementing them in any iterative search for the solution to the integral equation. As a consequence, it is crucial to have an interpolation routine for evaluating the function to be integrated at the irrational abscissas needed. This routine must be such as not to spoil the advantages of these integration methods. In this article we describe an interpolation routine that employs the Newtonforward, Newton-backward and Everett interpolation formulae. These formulae are displayed in section 2. In section 3 we describe the routine. In section 4 some general remarks are included.
.~
‘
f(x, +ph)
~
+
(P) k
~
k=
O~p~ 1,
~kj
(2)
the Newton-backward formula is
—
I p+k
~
fc~x +ph)—f+
~ i k= I
—
1
/i vkf
k
~‘
I ~p ~ 0
(3)
and the Everett interpolation formula is (n even) f(x
h)
+ s
+
P
k—I
f/I —p + k \ 2k + I
—
+
2k
)~
(~2kp ++k1
2kf
6
J
0~
~ I P
‘
(4) In these eqs. J~=f(x~); the combinatorial factors are defined by (P) = p (p I)... (p k + 1)/k! k and the finite differences are —
(5)
—
I.
forward
~f
=
~ (~ 1)~(k\f~~
~J/
j0 k
backward
vkf
central
(6a)
~‘
(~)
~ (— 1)~ f— 1=o 12k
,
(6b)
fs+k—r
(6c)
2k 82kf
~ (—l)~~ ~ j=0
J.A. Hernando
/ Newton
—
Everett interpolation of continuous functions
It must be noticed that the forward difference of order k contains all the function values between f, and fs+k’ the backward difference of order k contains the values between f~k and f~,and the central difference of order 2k, centred in J, contains the values between £—k and J.+k. Therefore, it is clear that the Everett interpolation is applicable only in the middle part of the interval and the Newton formulae are specially suited for the ends of the interval, The performance of these interpolation formulae is strongly dependent on the function to be interpolated, on the size of h and on the accuracy needed. This will be discussed in the next section. However, an important point must be emphasized: the tabular differences are a sensitive detector of roundoff and/or tabular errors [1,3]. In this discussion we consider only roundoff errors. When the regular flow of functional values is perturbed by random errors in some decimal place, it manifests itself in the appearance of an oscillating pattern in some tabular difference. The order of this finite difference tells us the accuracy of the functional values and the greatest finite difference that it makes sense to take into account.
3. The program The program is a double-precision function named FINT(XINT). The routine employs the Everett interpolation formula whenever possible and Newton formulae at the interval ends. If we look at the formulae (2), (3) and (4) we see that (for fixed n) we can rearrange them in a power series in p. In such a case each summand is a product of two factors: one solely dependent on the differences (i.e. on the subdivision of the interval and on the order of the interpolation) and the other being some power of p. Therefore, the factor dependent on the differences is evaluated only once in the first call to FINT. In all the subsequent calls these quantities are available. The distinction between these two cases is made with the logical command DIFF. When DIFF is true (first call), the quantities dependent on the differences are evaluated and control is returned to the main program without calculating any interpolated
75
value. In this case XINT is a dummy variable. When DIFF is false, the differences are supposed known and the program jumps to the direct evaluation of p and of the interpolated function at this point. It must be noticed that this methodology is efficient if we need to evaluate several interpolated values (whose number is at least equal to the number of known functional values NDIM, see below). In any other case it is more time consuming, than the more conventional method that evaluates exclusively the differences needed for each point and sums the series as it stands in (2), (3) or (4). There are other five logical commands: TWO, FOUR, SIX, EIGHT, TEN which direct the program to evaluate the function to the requested finite difference order. The last logical variable, WARN, is true when an error condition has been detected; in such a case the user must decide which action should be taken by the main program. All these logical commands are in the labelled common INT3. The functional values and the abscissas to which they correspond are in the vectors XCALC, FCALC. These two vectors and the separation between abscissas (H) are in the labelled common INTl. The dimension of these vectors (NDIM) and the integer LP are in the labelled common INT2. LP is an integer supplied by the user to indicate the output device where error messages should be sent. The error messages are of three kinds. The first one appears in the case that the five logical cornmands TWO,.. .,TEN are all false. The second message is for the case that the number of functional values (NDIM) does not suffice to evaluate the interpolation to the requested finite differences (NDIM must be greater than 2, 4, 6, 8, 10 for second, fourth, sixth, eighth and tenth interpolation order, respectively). The third error message is for the case that the point XINT, at which we want to evaluate the interpolated function, lies outside the interval covered by XCALC. There are several points worthy of mention with respect to error conditions. A first point concerns the change of interpolation orders (e.g. from TWO to FOUR). In this case all the quantities dependent on the differences should be recalculated;
76
J.A. Hernando
/ Newton —Everett interpolation of continuous functions
Table I CPU times and relative errors for the various options CPU time needed when: (for the Units, see text)
Relative error in the interpolated function when: _________________________________________________________________
f=exp(—x) cos(2~rx) DIFF=.T.
(T
1)
TWO FOUR SIX EIGHT TEN
II 24 48 78 120
___________________
___________________
(T,)
H=0.03
11=0.08
11=0.08
11=0.3
38 ‘U 43 46 48
tO ~ I0~ tO 8
10 iO~
to
to
lO~
~
IOb0
tO~
tO ~
tO
2
10 6 0 ~ 10 ~
10
~o
to
‘5
therefore, DIFF must be set true. A second point is that the routine does not check whether the vector XCALC is an ordered and equispaced vector; it assumes that to be true. If this is not true, the interpolation formulae (2), (3) and (4) do not apply. Another point concerns the possibility of having two commands true; e.g. FOUR and EIGHT. The routine does not detect this condition and performs the calculations ordered by the greatest difference command (EIGHT in this case). This superposition can be a source of errors, only some of which are covered by the checks included (e.g. that NDIM does not suffice for eighth order interpolation). With respect to the LP parameter: there is no default value, it must always be supplied by the user, The program is rather long and expensive in storage because our primary aim was to obtain an accurate and fast routine at the cost of increasing the storage requirements. Table 1 gives the CPU times for the various options that the program has built in. They are measured in units of the time needed to do a double precision multiplication, so that it can be estimated for another computer. If T1 is the time spent when DIFF is true (first call) and T2 when DIFF is false; the average time spent by interpolation in N calls to FINT is T
f=exp(—x)
DIFF=.F.
(NDIM*T 1 +N*T2 )/N.
In our system the time needed for a double precision exponential is 25 in these units, so it is seen that this routine is rather fast. In table 1 we also show the relative accuracy of
‘
to
the interpolation for two types of functions with two values of H. It is seen that we need relatively few values of FCALC to obtain a satisfactory accuracy.
4. Conclusions Very recently [4] a program has been published that evaluates the Tchebysheff series approximation to continuous functions. Tchebysheff series have a very attractive feature: they minimize the maximum error of the approximation (minimax approximation). Also interesting is the fact that to evaluate the Tchebysheff series one does not need an evenly spaced set of functional values. On the other hand, their coefficients must be determined by numerical integration and this implies both computation time and the dilemma of which integration method must be used. With the trapezoidal methods we lose accuracy and time and with the other methods an interpolation formalism in the spirit of this paper is needed. Therefore, it seems to us that these two routines can complement one another. To summarize, the main purpose of this paper was to present a fast and accurate interpolation routine. With an behaviour adequate choice of H (strongly dependent on the of the function, as can be seen in table I), it is possible to interpolate at practically the full accuracy of the given function. The economy in computation time increases with the number of interpolations calculated.
J.A. Hernando / Newton — Everett interpolation of continuous functions
References [I] A. Ralston, A first course in numerical analysis (McGrawHill, New York, 1965). [21 J.P. Hansen and I.R. McDonald, Theory of simple liquids (Academic Press, New York, 1976).
77
[3] Z. Kopal, Numerical analysis (Chapman & Hall, London, 1961). [41M.A. Christie and K.J.M. Moriarty, Comput. Phys. Commun. 23 (1981) 145.
78
J.A. Hernando
/
Newton — Everett interpolation of continuous functions
TEST RUN OUTPUT NEWTON—EVERETT INTERPOLATION FORMULAE FROM SECOND TO TENTII TABULAR DIFFLRENCF.S ORDER FUNCTION TO BE INTERPOLATED, FEXP(—X)COS(2PIX) 01FF ORDER TWO TWO TWO TWO
DIFFERENCE BETWEEN X VALUES, 9= 3.000—02
TWO
X VALUE 2.30000—02 l.4600D—0l 2.69i0D—01 3.92000—01 5.l5OiD—Ol 6.3800D—01 7.61000—01 8.8400D—01
FUNCTION VALUE= 9.670756160—01 5.25347650D—01 —9.100729170—02 —5.260101800—il —5.948488660—01 —3.418707650—01 3.226477460—02 3.001684980—il
INTERPOLATION= 9.669253700—01 5.25346730n—0l —9.100553800—02 —5.26005449D—01 —5.948419290—il —3. :18702900—il 3.2~560014D—02 3.081127220—01
RELATIVE ERROR= —l.554D—04 —1.7510—06 —1.9260—05 —8.9930—06 —1.1660—05 —1.3800—06 —2.719D—04 —1.810D—04
FOUR FOUR FOUR FOUR FOUR FOUR FOUR FOUR
2.30000—02 1.46000—01 2.6900D—01 3.92i0D—il 5.1500D—il 6.38000—01 7.61000—01 0.8400D—01
9.670756160—il 5.253476500—01 —9.10072917D—02 —5.260101800—01 —5.940488660—01 —3.41070765D—01 3.226477460—02 3.081684980—01
9.670784790—01 5.253476700—il —9.10072746D—02 —5.260101470—01 —5.948488330—01 —3.418707900—01 3.226468810—02 3.08168978D—01
2.96112—06 3.0210—08 —1.073D—i7 —6.229D—08 —5.5740—08 7.5210—08 —2.6790—06 1.5570—06
SIX SIX SIX SIX SIX SIX SIX SIX
2.30000—02 1.46000—01 2.69000—01 3.92000—01 5.lSOiD—0l 6.38000—01 7.61000—01 8.84000—il
9.670756160—01 5.253476500—01 —9.100729170—02 —5.260101800—01 —5.948488660—01 —3.418707650—01 3.226477460—02 3.081684980—01
9.670755570—01 5.253476500—01 —9.100729150—02 —5.26010l8iD—0l —5.948488660—01 —3.418707650—01 3.226477380—02 3.081684960—01
—6.1060—08 6.7610—10 —1.6880—09 —4.0980—10 —1.5280—10 1.2110—09 —2.4120—08 —6.6200—09
FlOE EIGH EIOH EIGH EIGH EIGH EIGH 12109
2.30000—02 1.46000—01 2.69000—01 3.92000—01 5. l500D—01 6.38000—01 7.61000—01 8.84000—01
9.670756160—01 5.2534765iD—01 —9.100729170—02 —5.260101800—01 —5.948488660—01 —3.418707650—01 3.226477460—02 3.081684980—01
9.670756170—01 5.253476500—01 —9.100729170—02 —5.260101800—01 —5.948408660—01 —3.418707650—01 3.226477460—02 3.081684970—01
1.1590—09 8.015D—12 —1.4150—11 —2.3620—12 1.2010—12 1.391D—1l —2.0090—10 —4.9520—10
TEN TEN TEN TEN TEN TEN TEN TEN
2.30000—02 1.46000—01 2.69000—01 3.92000—01 5.15000—01 6.38000—01 7.61000—01 8.84000—01
9.670756160—0]. 5.253476500—il —9.100729170—02 —5.260101800—01 —5.948488660—01 —3.418707650—01 3.226477460—02 3.081684980—01
9.670756160—01 5.253476500—il —9.l0i72917D—02 —5.261101800—01 —5.948488660—01 —3.41870765D—il. 3.226477450—02 3.081684980—01
—1.6760—11 5.8910—12 —1.1120—13 —9.89412—15 2.935D—14 1.3780—13 —5.3890—li 3.095D—ll
TWO
TWO TWO