A BASIC Program for Exact Curve Fitting

It is often required to find a polynomial to fit some function for which we know only a few values; the program given here in figure 1 calculates the coefficients of a polynomial such that the fit is exact for those data points.

If we have p+1 data points, then a polynomial of order p (or greater) can be found that will pass through all of these points (although in some special cases an exact polynomial of order less than p can be found). The minimum order polynomial is unique and it is the coefficients of this polynomial that are found by the program.

The program uses a method based on the Lagrange interpolating polynomial. For each value of abscissa: xj ( j is an integer: 0 less than or equal to j less than or equal to p ) we construct a polynomial: Aj (x) of order p that has the following properties: firstly, it has unity value at xj and secondly it is zero for all xi ( 0 less than or equal to i less than or equal to p; i not equal to j ). Aj (x) is constructed in two steps; firstly, using program statements: 130 - 190, a polynomial Fj (x) is constructed:

Fj (x) = ( x - x0 )( x - x1 ).......( x - xj-1 )( x - xj+1 )......( x - xp )

this clearly has the second property; by evaluating Fj (x) at xj (lines 220 - 230) we can find Aj (x) that has both properties (line 200):

    Fj (x)
Aj (x)   =  
    Fj ( evaluated at xj )

Having calculated an Aj (x) for all xj , weighting the coefficients of each Aj (x) by the data value corresponding to xj  and summing the polynomials gives the required function; this is calculated and printed by statements 275 - 315.

The program given here assumes that the data are given for values of x = 0,1,2,3.... but appropriate scaling will allow fitting to any set of equally spaced data (modifications to statement 160 will allow nonuniformly spaced data to be fitted). Although written and tested on a Hewlett Packard 85 computer all of the statements should be suitable for other interpreters with little or (in most cases) no modification.

By way of example, a fifth order fit to the function exp(x) is calculated, for x = 0,1,2,3,4 and 5, giving the coefficients of the required polynomial:

H0 = 1

H1 = 2.74952

H2 = -3.30606

H3 = 3.03500

H4 = -.885002

H5 =.124822

The error (H0 + H1 x + H2 x2 + H3 x3 + H4 x4 + H5 x5 - exp(x)) is plotted in figure 2; as can be seen the fit is exact at x = 0,1,2,3,4 and 5 and elsewhere the interpolation is good.

Lawrence J. Mayes
Department of Electrical Engineering Science
Essex University

Figure 1: Program Listing

100 COM A(57,56), C(57), H(56), Z(56)
105 INTEGER I, J, N, N1, P
110 CLEAR
115 LINPUT "Order of fit:",P$ @ P = VAL(P$) @ IF P > 56 THEN BEEP @ GOTO 115
120 GOSUB 245
125 FOR J = 0 TO P
130 GOSUB 235
135 C(0) = 1
140 FOR N = 0 TO P
145 N1 = N + 1
150 IF N = J THEN N = N + 1
155 IF N > P THEN 190
160 B0 = -N
165 FOR I = 1 TO N1
170 A(I,J) = B0 * C(I) + C(I-1)
175 NEXT I
180 A(0,J) = B0 * C(0)
185 FOR I = 0 TO N1 @ C(I) = A(I,J) @ NEXT I
190 NEXT N
195 GOSUB 220
200 FOR I = 0 TO P @ A(I,J) = A(I,J)/C(0) @ NEXT I
205 NEXT J
210 GOSUB 275
215 BEEP 100,190 @ DISP "finished" @ STOP
220 FOR I = P TO 1 STEP -1
225 C(I-1) = C(I) * J+C(I-1) @ NEXT I
230 RETURN
235 FOR I = 0 TO P + 1 @ C(I) = 0 @ NEXT I
240 RETURN
245 PRINT @ PRINT "INPUT DATA" @ PRINT "==========" @ PRINT
250 FOR I = 0 TO P
255 DISP @ DISP "Y value for abscissa value:"; I @ INPUT Z(I)
260 PRINT "Z("; I; ")= "; Z(I)
265 NEXT I
270 RETURN
275 PRINT @ PRINT "POLYNOMIAL FIT" @ PRINT "==============" @ PRINT
280 FOR I = 0 TO P @ H(I) = 0 @ NEXT I
285 FOR J = 0 TO P
290 FOR I = 0 TO P
295 H(I) = Z(J) * A(I,J) + H(I)
300 NEXT I
305 NEXT J
310 FOR I = 0 TO P @ PRINT "H("; I; ")= "; H(I) @ NEXT I
315 RETURN
320 END

error


NOTES

  • This article was published in ELECTRONIC ENGINEERING, 57, No 698; FEBRUARY 1985 p 38.
  • I have modified the programme to work with QBASIC - you can download a listing here.

Last updated: 1 Oct 2005;    © Lawrence Mayes, 1985/2000-2005