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