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: x_j ( j is an integer: 0 j p ) we construct a polynomial: A_j (x) of order p that has the following properties: firstly, it has unity value at x_j and secondly it is zero for all x_i ( 0 i p; i j ). A_j (x) is constructed in two steps; firstly, using program statements: 130 - 190, a polynomial F_j (x) is constructed:

F_j (x) = ( x - x₀ )( x - x₁ ).......( x - x_j-1 )( x - x_j+1 )......( x - x_p )

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

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

Having calculated an A_j (x) for all x_j , weighting the coefficients of each A_j (x) by the data value corresponding to x_j  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:

H₀ = 1

H₁ = 2.74952

H₂ = -3.30606

H₃ = 3.03500

H₄ = -.885002

H₅ =.124822

The error (H₀ + H₁ x + H₂ x² + H₃ x³ + H₄ x⁴ + H₅ x⁵ - 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

---

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.