'********************************************************* '* Approximation of a discrete real function F(x) by * '* least squares * '* ----------------------------------------------------- * '* Ref.: "Methodes de calcul numerique, Tome 2 by Claude * '* Nowakowski, PS1 Edition, 1984" [4]. * '* ----------------------------------------------------- * '* SAMPLE RUN: * '* * '* Number of points : 11 * '* * '* Degree of polynomial: 3 * '* * '* Function to approximate: * '* X(1), Y(1) = 0,0 * '* X(2), Y(2) = 0.1,0.2 * '* X(3), Y(3) = 0.2,0.4 * '* X(4), Y(4) = 0.3,0.6 * '* X(5), Y(5) = 0.4,0.8 * '* X(6), Y(6) = 0.5,1 * '* X(7), Y(7) = 0.6,0.8 * '* X(8), Y(8) = 0.7,0.6 * '* X(9), Y(9) = 0.8,0.4 * '* X(10), Y(10) = 0.9,0.2 * '* X(11), Y(11) = 1,0 * '* * '* Polynomial approximation of degree 3 (11 points) * '* Coefficients of polynomial: * '* A( 0 ) = -0.069930 * '* A( 1 ) = 3.496503 * '* A( 2 ) = -3.496503 * '* A( 3 ) = -0.000000 * '* * '* Approximated function: * '* X Y * '* 0.000000 -0.069930 * '* 0.100000 0.244755 * '* 0.200000 0.489510 * '* 0.300000 0.664336 * '* 0.400000 0.769231 * '* 0.500000 0.804196 * '* 0.600000 0.769231 * '* 0.700000 0.664336 * '* 0.800000 0.489510 * '* 0.900000 0.244755 * '* 1.000000 -0.069930 * '* * '* Basic Version By J-P Moreau, Paris. * '* (www.jpmoreau.fr) * '********************************************************* 'Program Approx DEFINT I-N DEFDBL A-H, O-Z NSIZE = 25 DIM C(NSIZE, NSIZE) DIM A(NSIZE), B(NSIZE), X(NSIZE), Xc(NSIZE), Y(NSIZE), Yx(NSIZE) CLS PRINT INPUT " Number of points : ", n n = n - 1 PRINT INPUT " Degree of polynomial: ", m n1 = n + 1: m1 = m + 1: m2 = m + 2 PRINT PRINT " Function to approximate:" FOR i = 1 TO n1 PRINT " X("; i; "), Y("; i; ") = "; INPUT "", X(i), Y(i) NEXT i FOR k = 1 TO m2 Xc(k) = 0# FOR i = 1 TO n1 Xc(k) = Xc(k) + X(i) ^ k NEXT i NEXT k yc = 0# FOR i = 1 TO n1 yc = yc + Y(i) NEXT i FOR k = 1 TO m Yx(k) = 0# FOR i = 1 TO n1 Yx(k) = Yx(k) + Y(i) * X(i) ^ k NEXT i NEXT k FOR i = 1 TO m1 FOR j = 1 TO m1 ij = i + j - 2 IF i = 1 AND j = 1 THEN C(1, 1) = n1 ELSE C(i, j) = Xc(ij) END IF NEXT j NEXT i B(1) = yc FOR i = 2 TO m1 B(i) = Yx(i - 1) NEXT i FOR k = 1 TO m FOR i = k + 1 TO m1 B(i) = B(i) - C(i, k) / C(k, k) * B(k) FOR j = k + 1 TO m1 C(i, j) = C(i, j) - C(i, k) / C(k, k) * C(k, j) NEXT j NEXT i NEXT k A(m1) = B(m1) / C(m1, m1) FOR i = m TO 1 STEP -1 s = 0# FOR k = i + 1 TO m1 s = s + C(i, k) * A(k) NEXT k A(i) = (B(i) - s) / C(i, i) NEXT i F\$ = "####.######" PRINT " Polynomial approximation of degree "; m; " ("; n + 1; " points)" PRINT PRINT " Coefficients of polynomial:" FOR i = 1 TO m1 PRINT " A("; i - 1; ") = "; : PRINT USING F\$; A(i) NEXT i PRINT PRINT " Approximated function:" PRINT " X Y" FOR i = 1 TO n1 xx = X(i): p = 0# FOR k = 1 TO m1 p = p * xx + A(m1 + 1 - k) NEXT k PRINT USING F\$; xx; : PRINT USING F\$; p NEXT i PRINT PRINT END ' end of file approx.bas