Curve fitting demo in
ABC FORTH
ChingTang Tseng
2015/4/12 Hamilton , New Zealand
1. Abstract
This article demonstrates it is quite simple to translate a typical traditional BASIC program into ABC FORTH
system. Original program implemented out curve fitting by Lagrange interpolating
polynomial method. A test program, according to the value of the exponential function
to be the input data, could give out all the coefficients of a polynomial by
the specified order.
2. Reference
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_{0}
)( x  x_{1} ).......( x  x_{j1} )( 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):
F_{j }(x)


A_{j }(x)

=


F_{j }( evaluated at x_{j} )

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_{0} = 1
H_{1} = 2.74952
H_{2} = 3.30606
H_{3} = 3.03500
H_{4} = .885002
H_{5} =.124822
The error (H_{0} + H_{1} x
+ H_{2} x^{2} + H_{3} x^{3} + H_{4} x^{4}
+ H_{5} x^{5}  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.
Department of Electrical
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(I1)
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(I1) = C(I) * J+C(I1) @ 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
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(I1)
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(I1) = C(I) * J+C(I1) @ 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.
Last
updated: 1 Oct
2005 ; © Lawrence Mayes,
1985/20002005
3. Program translated to the ABC FORTH system
\ ABC FORTH demo code:
Lagrange interpolating polynomial curve fitting.
\ Author: ChingTang Tseng,
Hamilton, New Zealand
\ Date：20120806(original),
20150409(renewed for posting)
\ Contact:
ilikeforth@gmail.com
\ Reference:
\ Lawrence J. Mayes
ELECTRONIC ENGINEERING, 57, No 698; FEBRUARY 1985 p 38.
\ Description:
\ This program uses a method
based on the Lagrange interpolating polynomial.
\ ZZ(i) is the input raw data
array.
\ Matrix AAA(i,j) is the
difference table manipulation matrix.
\ HH(i) is the output
coefficients array of the fitted polynomial.
\ Remark: Original array and
matrix structures in ABC FORTH system
\ without 0 indexes function.
0 index function to be included
\ after this program could be
implemented out.
5 INTEGERS I J N N1 P
: INPUTORDER BASIC
10 PRINT " Order of fit:
"
20 INPUTI P
30 IF P > 56 THEN 50
40 GOTO 70
50 PRINT " Order P cannot
> 56, Reenter please. "
60 GOTO 10
70 END ;
56 ARRAY ZZ
: SUB245 BASIC
10 RUN INPUTORDER
20 PRINT " Input data:
"
30 PRINT " ===========
"
40 FOR I = 0 TO P \ Index start from 0
50 PRINT " Y value for X
abscissa value: " ; I
60 INPUTR ZZ ( I )
70 PRINT " Y( " ; I
; " )=" ; { ZZ ( I ) }
80 NEXT I
90 END ;
57 ARRAY CC
: SUB235 BASIC
10 FOR I = 0 TO P + 1 \ Index is an arithmetical expression
20 LET { CC ( I ) = 0 }
30 NEXT I
40 END ;
56 ARRAY HH
57 56 MATRIX AAA
: SUB275 BASIC
10 PRINT " Polynomial
fit: "
20 PRINT "
=============== "
30 FOR I = 0 TO P
40 LET { HH ( I ) = 0 }
50 NEXT I
60 FOR J = 0 TO P
70 FOR I = 0 TO P
80 LET { HH ( I ) = ZZ ( J )
* AAA ( I J ) + HH ( I ) }
90 NEXT I
100 NEXT J
110 FOR I = 0 TO P
120 PRINT " HH( " ;
I ; " )= " ; { HH ( I ) }
130 NEXT I
140 END ;
INTEGER I1
: SUB220 BASIC
10 FOR I = P TO 1 STEP 1
20 LET I1 = I  1
30 LET { CC ( I1 ) = CC ( I
) * I>R ( J ) + CC ( I1 ) }
40 NEXT I
50 END ;
INTEGER Q
REAL B0
: MAIN BASIC
10 RUN SUB245
20 FOR J = 0 TO P
30 RUN SUB235
40 LET { CC ( 0 ) = 1 }
50 FOR N = 0 TO P
60 LET N1 = N + 1
70 IF N = J THEN 90
80 GOTO 100
90 LET N = N + 1
100 IF N > P THEN 200
105 LET Q = NEGATE N
110 LET { B0 = I>R ( Q ) }
120 FOR I = 1 TO N1
130 LET I1 = I  1
140 LET { AAA ( I J ) = B0 *
CC ( I ) + CC ( I1 ) }
150 NEXT I
160 LET { AAA ( 0 J ) = B0 *
CC ( 0 ) }
170 FOR I = 0 TO N1
180 LET { CC ( I ) = AAA ( I
J ) }
190 NEXT I
200 NEXT N
210 RUN SUB220
220 FOR I = 0 TO P
230 LET { AAA ( I J ) = AAA (
I J ) / CC ( 0 ) }
240 NEXT I
250 NEXT J
260 RUN SUB275
270 PRINT "
Finished."
280 END ;
\ Usage：
\ (1) Execute the following
instruction GIVEN first to get the input data.
\ (2) Next, execute MAIN,
then enter the data which you got above and
\ follows the Q & A procedures.
: GIVEN BASIC
10 RUN CR CR
20 FOR I = 0 TO 5
30 PRINT " Y( " ; I
; " )= " ; { EXP ( I>R ( I ) ) }
40 NEXT I
50 RUN CR CR
60 END ;
\S
Execution results:
GIVEN
Y( 0 )= 1.00000000000
Y( 1 )= 2.71828182846
Y( 2 )= 7.38905609893
Y( 3 )= 20.0855369232
Y( 4 )= 54.5981500331
Y( 5 )= 148.413159103
ok
MAIN
Order of fit: ? 5
Input data:
===========
Y value for X abscissa value:
0 ? 1.00000000000
Y( 0 )= 1.00000000000
Y value for X abscissa value:
1 ? 2.71828182846
Y( 1 )= 2.71828182846
Y value for X abscissa value:
2 ? 7.38905609893
Y( 2 )= 7.38905609893
Y value for X abscissa value:
3 ? 20.0855369232
Y( 3 )= 20.0855369232
Y value for X abscissa value:
4 ? 54.5981500331
Y( 4 )= 54.5981500331
Y value for X abscissa value:
5 ? 148.413159103
Y( 5 )= 148.413159103
Polynomial fit:
===============
HH( 0 )= 1.00000000000
HH( 1 )= 2.74952933754
HH( 2 )= 3.30606647675
HH( 3 )= 3.03499879102
HH( 4 )= .885001709373
HH( 5 )= .124821886021
Finished. ok
**************************************************
Original paper supplied
output result data are:
H0 = 1
H1 = 2.74952
H2 = 3.30606
H3 = 3.03500
H4 = .885002
H5 = .124822
**************************************************
4. This diagram shows the
error condition in the specified range.
5. All implementation are
finished in Win32Forth V6.14.03 addedin ABC FORTH V659 under W7 operating
system.