home *** CD-ROM | disk | FTP | other *** search
/ Programmer 7500 / MAX_PROGRAMMERS.iso / PASCAL / CURVFIT.ZIP / CURVFIT.PAS
Encoding:
Pascal/Delphi Source File  |  1993-01-04  |  4.0 KB  |  171 lines
  1. PROGRAM CURVFIT;
  2.  
  3. { The coefficients are found by minimizing the sum of the square
  4.   of the difference between the given data and the prediction.
  5.   The equation AG=B is solved by matrix inversion. An augmented
  6.   matrix A is created for matrix inversion by attaching B to the
  7.   last column of A. The Gaussian elimination method reduces the
  8.   original A to a unit matrix as the solution.
  9.  
  10.   This program was originally written in BASIC and presented in
  11.   MACHINE DESIGN Sept 6, 1984. It was rewritten in PASCAL by
  12.   Tim Kelly Jan 26, 1986.
  13. }
  14. Var
  15.   A          : array[0..25,0..26] of real;
  16.   C          : array[0..25,0..26] of real;
  17.   X,Y,W,Z    : array[0..25] of real;
  18.   N          : integer; { number of data points}
  19.   L          : integer; { order of polynomial}
  20.   M,R,I,J,K  : integer;
  21.   S          : real;
  22.   Dum        : Char;
  23.  
  24. {.........................................................................}
  25.  
  26. Procedure MATRIX;
  27.  
  28. BEGIN
  29.   For K:=1 to M do
  30.   Begin
  31.     S:=A[K,K];
  32.     For I:=K to R do A[K,I]:=A[K,I]/S;
  33.     For J:=1 to M do
  34.     If J <> K then
  35.     begin
  36.       S:=A[J,K];
  37.       For I:=K to R do A[J,I]:=A[J,I]-A[K,I]*S;
  38.     end;
  39.   End;
  40. END;
  41.  
  42. {.........................................................................}
  43.  
  44. Procedure ON(y,x: Integer);
  45. BEGIN
  46.   GotoXY(x,y);
  47.   Write('[]',^H);
  48. END;
  49.  
  50. {.........................................................................}
  51.  
  52. Procedure ONL(y,x: Integer);
  53. BEGIN
  54.   GotoXY(x,y);
  55.   Write('--',^H);
  56. END;
  57.  
  58. {.........................................................................}
  59.  
  60. Procedure OFF(x,y: Integer);
  61. BEGIN
  62.   GotoXY(x,y);
  63.   Write(' ',^H);
  64. END;
  65.  
  66. {.........................................................................}
  67.  
  68. Procedure BOX;
  69. Var
  70.   i,j : Integer;
  71. BEGIN
  72.   For i:=1 to 39 do
  73.   begin
  74.     j:=i*2;
  75.     ON(1,j); ONL(7,j);ON(24,j);
  76.   end;
  77.   For i:=1 to 24 do
  78.   begin
  79.     ON(i,2); ON(i,78);
  80.   end;
  81. END;
  82.  
  83. {.........................................................................}
  84.  Procedure OUTPUTCK;
  85.  BEGIN
  86.    Write(X[i]:10,' ',y[i]:10,' ',w[i]:10,' ',z[i]:10,' ');
  87.  END;
  88.  
  89. {=========================================================================}
  90.  
  91. BEGIN
  92.   clrscr; BOX;
  93.   gotoXY(30,3);write('GENERAL ORDER POLYNOMIAL');
  94.   gotoXY(33,4);write('CURVE FIT PROGRAM');
  95.   gotoXY(8,6);write('Number of data points :? ');Read(N);
  96.   gotoXY(47,6);write('Order of Polynomial :? ');Read(L);
  97.  
  98.   M:=L+1;  R:=L+2;
  99.  
  100.   For I:=1 to N do
  101.   Begin
  102.     gotoXY(35,12);Write('X(',I,')=? ');read(X[i]);
  103.     gotoXY(35,13);Write('Y(',I,')=? ');read(Y[i]);
  104.     For J:=35 to 60 do OFF(j,12);
  105.     For J:=35 to 60 do OFF(j,13);
  106.   End;
  107.  
  108.   For J:=1 to N do C[J,1]:=1;
  109.  
  110.   For I:=2 to M do
  111.   Begin
  112.     For J:=1 to N do C[J,I]:=C[J,I-1]*X[J];
  113.   End;
  114.  
  115.   For I:= 1 to M do
  116.   Begin
  117.     For J:=1 to M do
  118.     begin
  119.       A[I,J]:=0.;
  120.       For K:=1 to N do A[I,J]:=A[I,J]+C[K,I]*C[K,J];
  121.     end;
  122.   End;
  123.  
  124.   For I:=1 to M do
  125.   Begin
  126.     A[I,R]:=0;
  127.     For K:=1 to N do A[I,R]:=A[I,R]+C[K,I]*Y[K];
  128.   End;
  129.  
  130.   Matrix;
  131.  
  132.   For I:=1 to N do
  133.   Begin
  134.     W[I]:=0;
  135.     For J:=1 to M do
  136.     begin      { evaluate polynomial }
  137.       If J=1 then S:=1 else
  138.       begin
  139.         S:=X[I];
  140.         If J>2 then For K:=3 to J do S:=S* X[I];
  141.       end;
  142.       W[I]:=W[I]+A[J,R]* S;
  143.     end;
  144.     Z[I]:=Abs( W[I]-Y[I] ) / Y[I] *100.;
  145.   End;
  146.  
  147.   gotoXY(5,8);Write(' Coefficients :');
  148.   gotoXY(5,10);For I:=0 to L do write('    X^',I:1,'    ');
  149.   gotoXY(5,11);For I:=1 to M do Write(A[I,R]:10,' ');
  150.  
  151.   gotoXY(17,13);Write('     X          Y       Y(calc)   %(error)');
  152.   gotoXY(17,14);For I:=1 to 44 do Write('-');
  153.   For I:=1 to N do
  154.   begin
  155.     If I < 9
  156.     then begin
  157.       gotoXY(17,14+I); OUTPUTCK end
  158.     else begin
  159.       If I = 9
  160.       then begin
  161.          gotoXY(17,23);Write('Hit any return key to continue');Read(Dum);
  162.          For J:=17 to 70 do For K:=15 to 23 do OFF(J,K) ;
  163.          gotoXY(17,6+I);OUTPUTCK end
  164.       else begin
  165.       gotoXY(17,6+I);OUTPUTCK; end;
  166.     end;
  167.   end;
  168.   GotoXY(4,22);Read(dum);ClrScr;
  169. END.
  170.  
  171.