home *** CD-ROM | disk | FTP | other *** search
/ Programmer 7500 / MAX_PROGRAMMERS.iso / INFO / TURBOPAS / PAS_ENG.ZIP / LEAST6.PAS < prev    next >
Encoding:
Pascal/Delphi Source File  |  1985-07-19  |  3.9 KB  |  148 lines
  1. program least3;  { --> 226 }
  2. { Pascal program to perform a linear least-squares fit }
  3. { on the properties of steam with Gauss-Jordan routine }
  4. { Seperate modules needed:
  5.     GAUSSJ}
  6.  
  7.  
  8. const maxr = 20;  { data prints }
  9.  maxc = 4;  { polynomial terms }
  10.  
  11. type ary = array[1..maxr] of real;
  12.  arys = array[1..maxc] of real;
  13.  ary2 = array[1..maxr,1..maxc] of real;
  14.  ary2s = array[1..maxc,1..maxc] of real;
  15.  
  16. var p,t,v,
  17.  y,y_calc : ary;
  18.  resid  : ary;
  19.  coef,sig : arys;
  20.  nrow,ncol : integer;
  21.  correl_coef : real;
  22.  
  23.  
  24. procedure get_data(var p,t: ary; { independant variable }
  25.      var v:   ary; { dependant variable }
  26.      var nrow: integer); { length of vectors }
  27. { get values for n and arrays x,y }
  28.  
  29. var i : integer;
  30.  
  31. begin
  32.   nrow:=12;
  33.   t[1]:=400; p[1]:=120; v[1]:=4.079;
  34.   t[2]:=450; p[2]:=120; v[2]:=4.36;
  35.   t[3]:=500; p[3]:=120; v[3]:=4.633;
  36.   t[4]:=400; p[4]:=140; v[4]:=3.466;
  37.   t[5]:=450; p[5]:=140; v[5]:=3.713;
  38.   t[6]:=500; p[6]:=140; v[6]:=3.952;
  39.   t[7]:=400; p[7]:=160; v[7]:=3.007;
  40.   t[8]:=450; p[8]:=160; v[8]:=3.228;
  41.   t[9]:=500; p[9]:=160; v[9]:=3.440;
  42.   t[10]:=400; p[10]:=180; v[10]:=2.648;
  43.   t[11]:=450; p[11]:=180; v[11]:=2.850;
  44.   t[12]:=500; p[12]:=180; v[12]:=3.042;
  45.   for i:=1 to nrow do
  46.     t[i]:=t[i]+460.0 { convert to Rankine }
  47. end;  { proceddure get data }
  48.  
  49. procedure write_data;
  50. { print out the answers }
  51. var i : integer;
  52. begin
  53.   writeln;
  54.   writeln('  I      P      T      V        Y    YCALC     %RES');
  55.   for i:=1 to nrow do
  56.     writeln(i:3,p[i]:7:1,t[i]:7:1,v[i]:7:3,y[i]:9:2,y_calc[i]:9:2,
  57.   (100.0*resid[i]/y[i]):9:2);
  58.   writeln; writeln(' Coefficients errors ');
  59.   writeln(coef[1],' ',sig[1],' Constant term');
  60.   for i:=2 to ncol do
  61.     writeln(coef[i],' ',sig[i]);  { other terms }
  62.   writeln;
  63.   writeln('Correlation coefficient is ',correl_coef:8:5)
  64. end;  { write_data }
  65.  
  66. {procedure square(x: ary2;
  67.    y: ary;
  68.       var a: ary2s;
  69.       var g: arys;
  70.   nrow,ncol: integer);}
  71. { matrix multiplication routine }
  72. { a= transpose x times x }
  73. { g= y times x }
  74. {$I SQUARE.LIB }
  75.  
  76. {external procedure gaussj(var b: ary2s;
  77.          y: arys;
  78.      var coef: arys;
  79.      ncol:  integer;
  80.      var error: boolean);
  81. }
  82. {$I GAUSSJ.LIB }
  83.  
  84. procedure linfit(p,t,v: ary; { independant variable }
  85.    var y: ary; { dependent variable }
  86.    var y_calc: ary; { calculated dep. variable }
  87.    var resid:  ary; { array of residuals }
  88.    var coef:   arys; { coefficients }
  89.    var sig:    arys; { error on coefficients }
  90.    nrow:       integer; { length of array }
  91.    var ncol:   integer); { number of terms }
  92.  
  93. { least squares fit to nrow sets of x and y pairs of points }
  94. { Seperate procedures needed:
  95.  SQUARE -> form square coefficient matrix
  96.  GAUSSJ -> Gauss-Jordan elimination }
  97.  
  98. const r = 85.76; { gas constant for steam }
  99.  
  100. var xmatr  : ary2;  { data matrix }
  101.  a  : ary2s; { coefficient matrix }
  102.  g  : arys;  { constant vector }
  103.  error  : boolean;
  104.  i,j,nm  : integer;
  105.  power,yi,yc,srs,see,
  106.  sum_y,sum_y2 : real;
  107.  
  108. begin  { procedure linfit }
  109.   ncol:=2; { number of terms }
  110.   for i:=1 to nrow do
  111.     begin  { setup matrix }
  112.       power:=t[i];
  113.       xmatr[i,1]:=p[i]/power; { first column }
  114.       xmatr[i,2]:=sqrt(p[i]); { second column }
  115.       y[i]:=v[i]*p[i]-r*t[i]/144.0
  116.     end;
  117.   square(xmatr,y,a,g,nrow,ncol);
  118.   gaussj(a,g,coef,ncol,error);
  119.   sum_y:=0.0;
  120.   sum_y2:=0.0;
  121.   srs:=0.0;
  122.   for i:=1 to nrow do
  123.     begin
  124.       yi:=y[i];
  125.       yc:=0.0;
  126.       for j:=1 to ncol do
  127.  yc:=yc+coef[j]*xmatr[i,j];
  128.       y_calc[i]:=yc;
  129.       resid[i]:=yc-yi;
  130.       srs:=srs+sqr(resid[i]);
  131.       sum_y:=sum_y+yi;
  132.       sum_y2:=sum_y2+yi*yi
  133.     end;
  134.   correl_coef:=sqrt(1.0-srs/(sum_y2-sqr(sum_y)/nrow));
  135.   if nrow=ncol then nm:=1
  136.   else nm:=nrow-ncol;
  137.   see:=sqrt(srs/nm);
  138.   for i:=1 to ncol do  { errors on solution }
  139.     sig[i]:=see*sqrt(a[i,i])
  140. end; { linfit }
  141.  
  142. begin  { main program }
  143.   clrscr;
  144.   get_data(p,t,v,nrow);
  145.   linfit(p,t,v,y,y_calc,resid,coef,sig,nrow,ncol);
  146.   write_data
  147. end.
  148.