Category : Pascal Source Code
Archive   : PAS-SCI.ZIP
Filename : NLIN3.PAS
{ Pascal program to perform a nonlinear least-squares fit for the diffusion
of Zn in CU }
const maxr = 20; { data prints }
maxc = 4; { polynomial terms }
r = 1.987; { gas constant }
type
index = 1..maxr;
ary = array[index] of real;
arys = array[1..maxc] of real;
ary2 = array[1..maxr,1..maxc] of real;
var
x,y,y_calc : ary;
t,d,ex : ary;
coef : arys;
i,n : integer;
nrow,ncol : integer;
done,error : boolean;
correl_coef,srs,
a,b,x2 : real;
procedure get_data(var x,y: ary;
var n: integer);
{ get values for n and arrays t,d }
var i : integer;
begin
n:=7;
t[1]:=600.0; d[1]:=1.4E-12;
t[2]:=650.0; d[2]:=5.5E-12;
t[3]:=700.0; d[3]:=1.8E-11;
t[4]:=750.0; d[4]:=6.1E-11;
t[5]:=800.0; d[5]:=1.6E-10;
t[6]:=850.0; d[6]:=4.4E-10;
t[7]:=900.0; d[7]:=1.2E-9;
for i:=1 to n do
begin
x[i]:=1.0/(t[i]+273.0);
y[i]:=d[i]
end
end; { proceddure get data }
procedure write_data;
{ print out the answers }
var i : integer;
begin
writeln;
writeln;
writeln(' I TC D DCALC');
for i:=1 to n do
writeln(i:3,t[i]:8:0,d[i],' ',y_calc[i]);
writeln; writeln(' Coefficients ');
writeln(coef[1],' constant term');
for i:=2 to ncol do
writeln(coef[i]); { other terms }
writeln;
writeln('D0=',a:7:2,' cm sq/sec.');
writeln('Q =',(-r*b/1000.0):8:2,'kcal/mole');
writeln;writeln('SRS= ',srs:8:4)
end; { write_data }
procedure func(b: real;
var fb,dfb: real);
var i : integer;
s1,s2,s3,s4,s5,s6,
ex1,ex2,xi,
x2,yi,y2 : real;
begin
s1:=0.0;
s2:=0.0;
s3:=0.0;
s4:=0.0;
s5:=0.0;
s6:=0.0;
for i:=1 to n do
begin
xi:=x[i];
x2:=xi*xi;
yi:=y[i];
y2:=yi*yi;
ex1:=exp(b*xi);
ex[i]:=ex1;
ex2:=ex1*ex1;
s1:=s1+xi*ex2/y2;
s2:=s2+ex1/yi;
s3:=s3+xi*ex1/yi;
s4:=s4+ex2/y2;
s5:=s5+2.0*x2*ex2/y2;
s6:=s6+x2*ex1/yi
end;
fb:=s1*s2-s3*s4;
dfb:=s2*s5-s1*s3-s4*s6;
a:=s2/s4
end; { func }
procedure newton(var x: real);
const tol = 1.0E-6;
max = 20;
var fx,dfx,dx,x1 : real;
i : integer;
begin { newton }
error:=false;
i:=0;
repeat
i:=i+1;
x1:=x;
func(x,fx,dfx);
if dfx=0.0 then
begin
error:=true;
x:=1.0;
writeln('ERROR: slope zero')
end
else
begin
dx:=fx/dfx;
x:=x1-dx;
end
until
error or
(i>max) or
(abs(dx)<=abs(tol*x));
if i>max then
begin
writeln(chr(7),'ERROR: no convergence in ',max,' loops');
error:=true
end
end; { newton }
procedure nlin(x,y: ary;
var y_calc: ary;
n: integer);
{ fits the diffusion equation through n sets of x and y pairs of points }
var
resid : ary;
matr : ary2;
i : integer;
xi,yi,sum_x,
sum_y,sum_y2,b1,
sum_xy,sum_x2 : real;
begin { nlin }
ncol:=2; { number of terms }
sum_x:=0.0;
sum_y:=0.0;
sum_xy:=0.0;
sum_x2:=0.0;
for i:=1 to n do
begin
xi:=x[i];
yi:=ln(y[i]);
sum_x:=sum_x+xi;
sum_y:=sum_y+yi;
sum_y2:=sum_y2+yi*yi;
sum_xy:=sum_xy+xi*yi;
sum_x2:=sum_x2+xi*xi
end;
b:=(sum_xy-sum_x*sum_y/n)/(sum_x2-sqr(sum_x)/n);
newton(b);
coef[1]:=a;
coef[2]:=b;
srs:=0.0;
for i:=1 to n do
begin
y_calc[i]:=a*ex[i];
if y[i]<>0.0 then
resid[i]:=y_calc[i]/y[i]-1.0
else resid[i]:=y[i]/y_calc[i]-1.0;
srs:=srs+sqr(resid[i])
end
end; { nlin }
begin { main program }
ClrScr;
get_data(x,y,n);
nlin(x,y,y_calc,n);
write_data
end.
Very nice! Thank you for this wonderful archive. I wonder why I found it only now. Long live the BBS file archives!
This is so awesome! 😀 I’d be cool if you could download an entire archive of this at once, though.
But one thing that puzzles me is the “mtswslnkmcjklsdlsbdmMICROSOFT” string. There is an article about it here. It is definitely worth a read: http://www.os2museum.com/wp/mtswslnk/