home *** CD-ROM | disk | FTP | other *** search
/ Programmer 7500 / MAX_PROGRAMMERS.iso / INFO / TURBOPAS / PAS_ENG.ZIP / NLIN3.PAS < prev    next >
Encoding:
Pascal/Delphi Source File  |  1985-07-19  |  3.9 KB  |  192 lines
  1. { die ln-funktion bei MT+ ist fuer Werte kleiner 1E-5 aeussertst langsam,
  2.   so dass man das Gefuehl bekommt, das dass Programm sich aufhaengt. In-
  3.   folgedessen konnte ich dieses Programm nicht verifizieren. -Juergen }
  4.  
  5. program nlin3;  { -> 310 }
  6. { Pascal program to perform a nonlinear least-squares fit for the diffusion
  7.  of Zn in CU }
  8.  
  9. const maxr = 20;  { data prints }
  10.  maxc = 4;  { polynomial terms }
  11.  r = 1.987; { gas constant }
  12. type
  13.  index = 1..maxr;
  14.  ary = array[index] of real;
  15.  arys = array[1..maxc] of real;
  16.  ary2 = array[1..maxr,1..maxc] of real;
  17.  
  18. var
  19.  x,y,y_calc : ary;
  20.  t,d,ex  : ary;
  21.  coef  : arys;
  22.  i,n  : integer;
  23.  nrow,ncol : integer;
  24.  done,error : boolean;
  25.  correl_coef,srs,
  26.  a,b,x2  : real;
  27.  
  28.  
  29. procedure get_data(var x,y: ary;
  30.      var n: integer);
  31. { get values for n and arrays t,d }
  32.  
  33. var i : integer;
  34.  
  35. begin
  36.   n:=7;
  37.   t[1]:=600.0; d[1]:=1.4E-12;
  38.   t[2]:=650.0; d[2]:=5.5E-12;
  39.   t[3]:=700.0; d[3]:=1.8E-11;
  40.   t[4]:=750.0; d[4]:=6.1E-11;
  41.   t[5]:=800.0; d[5]:=1.6E-10;
  42.   t[6]:=850.0; d[6]:=4.4E-10;
  43.   t[7]:=900.0; d[7]:=1.2E-9;
  44.   for i:=1 to n do
  45.     begin
  46.       x[i]:=1.0/(t[i]+273.0);
  47.       y[i]:=d[i]
  48.     end
  49. end;  { proceddure get data }
  50.  
  51. procedure write_data;
  52. { print out the answers }
  53. var i : integer;
  54. begin
  55.   writeln;
  56.   writeln;
  57.   writeln('     I      TC     D      DCALC');
  58.   for i:=1 to n do
  59.     writeln(i:3,t[i]:8:0,d[i],' ',y_calc[i]);
  60.   writeln; writeln(' Coefficients ');
  61.   writeln(coef[1],' constant term');
  62.   for i:=2 to ncol do
  63.     writeln(coef[i]);  { other terms }
  64.   writeln;
  65.   writeln('D0=',a:7:2,' cm sq/sec.');
  66.   writeln('Q =',(-r*b/1000.0):8:2,'kcal/mole');
  67.   writeln;writeln('SRS= ',srs:8:4)
  68. end;  { write_data }
  69.  
  70. procedure func(b: real;
  71.       var fb,dfb: real);
  72. var i  : integer;
  73.  s1,s2,s3,s4,s5,s6,
  74.  ex1,ex2,xi,
  75.  x2,yi,y2 : real;
  76. begin
  77.   s1:=0.0;
  78.   s2:=0.0;
  79.   s3:=0.0;
  80.   s4:=0.0;
  81.   s5:=0.0;
  82.   s6:=0.0;
  83.     for i:=1 to n do
  84.       begin
  85.  xi:=x[i];
  86.  x2:=xi*xi;
  87.  yi:=y[i];
  88.  y2:=yi*yi;
  89.  ex1:=exp(b*xi);
  90.  ex[i]:=ex1;
  91.  ex2:=ex1*ex1;
  92.  s1:=s1+xi*ex2/y2;
  93.  s2:=s2+ex1/yi;
  94.  s3:=s3+xi*ex1/yi;
  95.  s4:=s4+ex2/y2;
  96.  s5:=s5+2.0*x2*ex2/y2;
  97.  s6:=s6+x2*ex1/yi
  98.       end;
  99.     fb:=s1*s2-s3*s4;
  100.     dfb:=s2*s5-s1*s3-s4*s6;
  101.     a:=s2/s4
  102. end;  { func }
  103.  
  104.  
  105. procedure newton(var x: real);
  106. const tol  = 1.0E-6;
  107.  max  = 20;
  108. var fx,dfx,dx,x1 : real;
  109.  i  : integer;
  110.  
  111. begin { newton }
  112.   error:=false;
  113.   i:=0;
  114.   repeat
  115.     i:=i+1;
  116.     x1:=x;
  117.     func(x,fx,dfx);
  118.     if dfx=0.0 then
  119.       begin
  120.  error:=true;
  121.  x:=1.0;
  122.  writeln('ERROR: slope zero')
  123.       end
  124.     else
  125.       begin
  126.  dx:=fx/dfx;
  127.  x:=x1-dx;
  128.       end
  129.   until
  130.     error or
  131.       (i>max) or
  132.  (abs(dx)<=abs(tol*x));
  133.   if i>max then
  134.     begin
  135.       writeln(chr(7),'ERROR: no convergence in ',max,' loops');
  136.       error:=true
  137.     end
  138. end;  { newton }
  139.  
  140. procedure nlin(x,y: ary;
  141.  var y_calc: ary;
  142.         n: integer);
  143. { fits the diffusion equation through n sets of x and y pairs of points }
  144. var
  145.  resid  : ary;
  146.  matr  : ary2;
  147.  i  : integer;
  148.  xi,yi,sum_x,
  149.  sum_y,sum_y2,b1,
  150.  sum_xy,sum_x2 : real;
  151. begin  { nlin }
  152.   ncol:=2; { number of terms }
  153.   sum_x:=0.0;
  154.   sum_y:=0.0;
  155.   sum_xy:=0.0;
  156.   sum_x2:=0.0;
  157.   for i:=1 to n do
  158.     begin
  159.       xi:=x[i];
  160.       yi:=ln(y[i]);
  161.       sum_x:=sum_x+xi;
  162.       sum_y:=sum_y+yi;
  163.       sum_y2:=sum_y2+yi*yi;
  164.       sum_xy:=sum_xy+xi*yi;
  165.       sum_x2:=sum_x2+xi*xi
  166.     end;
  167.   b:=(sum_xy-sum_x*sum_y/n)/(sum_x2-sqr(sum_x)/n);
  168.   newton(b);
  169.   coef[1]:=a;
  170.   coef[2]:=b;
  171.   srs:=0.0;
  172.   for i:=1 to n do
  173.     begin
  174.       y_calc[i]:=a*ex[i];
  175.       if y[i]<>0.0 then
  176.  resid[i]:=y_calc[i]/y[i]-1.0
  177.       else resid[i]:=y[i]/y_calc[i]-1.0;
  178.       srs:=srs+sqr(resid[i])
  179.     end
  180. end;   { nlin }
  181.  
  182.  
  183. begin  { main program }
  184.   clrscr;
  185. writeln(' start get_data ');
  186.   get_data(x,y,n);
  187. writeln(' start nlin ');
  188.   nlin(x,y,y_calc,n);
  189. writeln(' start write_data ');
  190.   write_data
  191. end.
  192.