## 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

data points, then a polynomial of orderp+1(or greater) can be found that will pass through all of these points (although in some special cases an exact polynomial of order less thanpcan be found). The minimum order polynomial is unique and it is the coefficients of this polynomial that are found by the program.pThe program uses a method based on the Lagrange interpolating polynomial. For each value of abscissa:

(x_{j}is an integer: 0jj) we construct a polynomial:pof orderA_{j }(x)that has the following properties: firstly, it has unity value atpand secondly it is zero for allx_{j}( 0x_{i}i;pi).jis constructed in two steps; firstly, using program statements: 130 - 190, a polynomialA_{j }(x)is constructed:F_{j }(x)

=F_{j }(x)( x - x_{0})( x - x_{1}).......( x - x_{j-1})( x - x_{j+1})......( x - x_{p})this clearly has the second property; by evaluating

atF_{j }(x)(lines 220 - 230) we can findx_{j}that has both properties (line 200):A_{j }(x)

F_{j }(x)A_{j }(x)= F_{j }( evaluated at x_{j})

Having calculated an

for allA_{j }(x), weighting the coefficients of eachx_{j}by the data value corresponding toA_{j }(x)and summing the polynomials gives the required function; this is calculated and printed by statements 275 - 315.x_{j}The program given here assumes that the data are given for values of

= 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.xBy way of example, a fifth order fit to the function

is calculated, forexp(x)= 0,1,2,3,4 and 5, giving the coefficients of the required polynomial:x

= 1H_{0}

= 2.74952H_{1}

= -3.30606H_{2}

= 3.03500H_{3}

= -.885002H_{4}

=.124822H_{5}The error (

-H_{0}+ H_{1}x + H_{2}x^{2}+ H_{3}x^{3}+ H_{4}x^{4}+ H_{5}x^{5}) is plotted in figure 2; as can be seen the fit is exact atexp(x)= 0,1,2,3,4 and 5 and elsewhere the interpolation is good.xLawrence 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.

Electronics Index |
Feedback Form |