程序如下。每次运行总是报：??? Error: File: E:\work\K4.m Line: 134 Column: 1
Illegal use of reserved keyword "function".
我开始以为是子程序出了问题，可是把每个子程序提出来单独编程却能通过，奇怪啊。请高人指点一二

function K4
tic
clc;format short;
n=4;
Pt=1;
R=83.145;
%%--------------------------------------------------------------------
A=[1 0 0 0 0 1.4311 1.432;
0 0 1 0 0 0.92 1.4;
0 1 0 1 0 0.9011 0.848;
0 1 0 2 0 0.6744 0.54;
0 1 0 1 0 1 1.2;
0 0 0 0 0 0 0;
0 0 0 0 0 0 0;
0 0 0 0 0 0 0;
0 0 0 0 0 0 0;
0 0 0 0 0 0 0;];
B=[0 -181 16.51 16.51 249.1 0 0 0 0 0;
289.6 0 300 300 -229.1 0 0 0 0 0;
697.2 1318 0 0 986.5 0 0 0 0 0;
697.2 1318 0 0 986.5 0 0 0 0 0;
-137.1 353.5 156.4 156.4 0 0 0 0 0 0;
0 0 0 0 0 0 0 0 0 0;
0 0 0 0 0 0 0 0 0 0;
0 0 0 0 0 0 0 0 0 0;
0 0 0 0 0 0 0 0 0 0;
0 0 0 0 0 0 0 0 0 0;];
D=sum(A);
mm=zeros(10,1);
for i=1:10
m=0;
for j=1:n
m=A(i,j)+m;
end
if m>0
mm(i)=100;
else
mm(i)=i;
end
end
N=min(mm)-1;
%--------------------------------------------------------------------------------------------
CriticalParameter=[8.094 512.6 118 0.224 0.559 0.0878 0.056;
6.382 516.2 167 0.248 0.635 0.0878 0.0572;
22.048 647.1 56 0.229 0.344 0.0279 0.0229;
5.175 536.78 219 0.254 0.629 0.0878 0.0439;
0 0 0 0 0 0 0;];
CriticalTemperature=[512.600 514.397 518.343 524.551 0.00;
0.000 516.200 520.160 526.389 0.00;
0.000 0.000 647.100 530.428 0.00;
0.000 0.000 0.000 536.780 0.00;
0.000 0.000 0.000 0.000 0.000;];
Antoine=[23.4804 23.8047 23.1964 22.4367 0.0000;
3626.5500 3803.9800 3816.4400 3166.3800 0.0000;
-34.2900 -41.6800 -46.1300 -80.1500 0.0000;];
%-------------------------------------------------------
K=zeros(1000,9);
xx=zeros(1,n);
num=0;
xx(1)=0.5;
xx(2)=0.2;
xx(3)=0.15;
xx(4)=0.15;
if sum(xx)==1
num=num+1;
t0=zeros(1,n);
for i=1:n
t0(i)=Antoine(2,i)/(Antoine(1,i)-lg(100000*Pt))-Antoine(3,i);
end
nn=1;
T=zeros(1,100);
T(nn)=Temper(xx,t0,n,Pt,R,Antoine,A,B,D,N);
Q=ones(1,n);
F=zeros(1,100);
y=zeros(1,n);
%-----------------------------------------------------------
while nn<50
yc=ActivityCoefficient(xx,T(nn),n,R,A,B,D,N);
VPsat=VaporPressure(T(nn),n,Antoine);
Qv=VP(T(nn),n,Pt,R,CriticalParameter,CriticalTemperature,Antoine);
for i=1:n
y(i)=xx(i)*yc(i)*VPsat(i)*Qv(i)/(Pt*Q(i));
end
SUMY0=sum(y);
M=1;
while M<50
Qm=mixVP(y,T(nn),n,Pt,R,CriticalParameter,CriticalTemperature,Antoine);
for i=1:n
y(i)=xx(i)*yc(i)*VPsat(i)*Qv(i)/(Pt*Qm(i));
end
SUMY=sum(y);
if abs((SUMY-SUMY0)/SUMY0)>0.00001
M=M+1;
SUMY0=SUMY;
y(i)=y(i)/SUMY;
else
break
end
end
F(nn)=SUMY-1;
%-----------------------------------------------------------
if abs(F(nn))>0.000001
nn=nn+1;
if nn<=2
if SUMY>1
T(nn)=0.98*T(nn-1);
else
T(nn)=1.02*T(nn-1);
end
else
dT=-(T(nn-1)-T(nn-2))*F(nn-1)/(F(nn-1)-F(nn-2));
T(nn)=T(nn-1)+dT;
if abs(dT)<0.01
break
end
end
else
break
end
%----------------------------------------------------------
K=zeros(1,n);
for i=1:n
y(i)=y(i)/SUMY;
K(i)=y(i)/x(i);
end
Tb=T(nn);
disp(y);
disp(K);
disp(Tb);
end
TIME=toc
disp(TIME);
%----------------------------------------------------------------------------
function f=coefficient(t,n,R,CriticalParameter,CriticalTemperature,Antoine)
Q=zeros(1,5);
Bij=zeros(5,5);
P=zeros(5,5);
W=zeros(5,5);
aij=zeros(5,5);
bij=zeros(5,5);
J0=zeros(5,5);
J1=zeros(5,5);
J2=zeros(5,5);
Tr=zeros(5,5);
XX=zeros(1,5);
for i=1:n
for j=i:n
Tr(i,j)=t/CriticalTemperature(i,j);
if i==j
P(i,j)=10*CriticalParameter(i,1);
W(i,j)=CriticalParameter(i,5);
aij(i,j)=CriticalParameter(i,6);
bij(i,j)=CriticalParameter(i,7);
else
P(i,j)=4*CriticalTemperature(i,j)*(10*CriticalParameter(i,1)*CriticalParameter(i,3)/CriticalParameter(i,2)+10*CriticalParameter(j,1)*CriticalParameter(j,3)/CriticalParameter(j,2))/(CriticalParameter(i,3)^(1/3)+CriticalParameter(j,3)^(1/3))^3;
W(i,j)=0.5*(CriticalParameter(i,5)+CriticalParameter(j,5));
aij(i,j)=0.5*(CriticalParameter(i,6)+CriticalParameter(j,6));
bij(i,j)=0.5*(CriticalParameter(i,7)+CriticalParameter(j,7));
end
J0(i,j)=0.1445-0.33/Tr(i,j)-0.1385/Tr(i,j)^2-0.0121/Tr(i,j)^3-0.000607/Tr(i,j)^8;
J1(i,j)=0.0637+0.331/Tr(i,j)^2-0.423/Tr(i,j)^3-0.008/Tr(i,j)^8;
J2(i,j)=aij(i,j)/Tr(i,j)^6-bij(i,j)/Tr(i,j)^8;
Bij(i,j)=(J0(i,j)+W(i,j)*J1(i,j)+J2(i,j))*R*CriticalTemperature(i,j)/P(i,j);
end
end
f=Bij;
%-----------------------------------------------------------------------------------------
function f=mixVP(y,t,n,Pt,R,CriticalParameter,CriticalTemperature,Antoine)
BB=0;
Bij=coefficient(t,n,R,CriticalParameter,CriticalTemperature,Antoine);
for i=1:n
for j=i:n
BB=BB+y(i)*y(j)*Bij(i,j);
end
end
yy=[Pt -R*t -BB*R*t];
Vx=roots(yy);
Z=1+BB/max(Vx);
for i=1:n
for j=1:n
Q(i)=Q(i)+y(j)*Bij(i,j);
end
Q(i)=exp(2/max(Vx)*Q(i)-log(Z));
end
f=Q;
%---------------------------------------------------------------------------------------------
function f=VaporPressure(t,n,Antoine)
Pv=zeros(1,n);
for i=1:n
Pv(i)=exp(Antoine(1,i)-Antoine(2,i)/(Antoine(3,i)+t))/100000;
end
f=Pv;
%--------------------------------------------------------------------------------------------------
function f=VP(t,n,Pt,R,CriticalParameter,CriticalTemperature,Antoine)
Bij=coefficient(t,n,R,CriticalParameter,CriticalTemperature,Antoine);
P=VaporPressure(t,n,Antoine);
Qv=zeros(1,n);
Zz=zeros(1,n);
Vvx=zeros(1,n);
for i=1:n
yyy=[P(i) -R*t -Bij(i,i)*R*t];
Vvx=roots(yyy);
Zz(i)=1+Bij(i,i)/max(Vvx);
Qv(i)=exp(2/max(Vvx)*Bij(i,i)-log(Zz(i)));
end
f=Qv;
%---------------------------------------------------------------------------------------------
function f=ActivityCoefficient(xx,t,n,R,A,B,D,N)
C=zeros(N,N);
for i=1:N
for j=1:N
if i==j
C(i,j)=1;
else
C(i,j)=exp(-1*B(i,j)/t);
end
end
end
r=zeros(1,n);
q=zeros(1,n);
v=zeros(1,n);
s=zeros(1,n);
l=zeros(1,n);
for j=1:n
r(j)=0;
q(j)=0;
for i=1:N
r(j)=r(j)+A(i,j)*A(i,6);
q(j)=q(j)+A(i,j)*A(i,7);
end
end
for j=1:n
v(j)=xx(j)*r(j)/dot(xx,r);
s(j)=xx(j)*q(j)/dot(xx,q);
l(j)=5*(r(j)-q(j))-r(j)+1;
end
lny2=zeros(n,1);
sum1=0;
for j=1:n
sum1=sum1+xx(j)*l(j);
end
for j=1:n
lny2(j)=log(v(j)/xx(j))+5*q(j)*log(s(j)/v(j))+l(j)-v(j)/xx(j)*sum1;
end
X=zeros(N,6);
for j=1:n
for i=1:N
X(i,j)=A(i,j)/D(j);
end
end
sumxD=0;
for j=1:n
sumxD=sumxD+xx(j)*D(j);
end
for i=1:N
sumxA=0;
for j=1:n
sumxA=sumxA+xx(j)*A(i,j);
end
X(i,6)=sumxA/sumxD;
end
SX=zeros(N,6);
for j=1:n
for i=1:N
SX(i,j)=X(i,j)*A(i,7)/dot(X(:,j),A(1:N,7));
end
end
for i=1:N
SX(i,6)=X(i,6)*A(i,7)/dot(X(:,6),A(1:N,7));
end
E=zeros(N,6);
for j=1:n
for i=1:N
if SX(i,j)>0
E(i,j)=dot(SX(:,j),C(:,i));
else
E(i,j)=0;
end
end
end
for i=1:N
E(i,6)=dot(SX(:,6),C(:,i));
end
F=zeros(N,6);
f=zeros(N,6);
ff=zeros(1,10);
for j=1:6;
for i=1:N;
if E(i,j)>0
f(i,j)=SX(i,j)/E(i,j);
else
f(i,j)=0;
end
end
end
for j=1:n;
for i=1:N;
if E(i,j)>0
F(i,j)=dot(f(:,j),C(i,);
else
F(i,j)=0;
end
end
end
for i=1:N
F(i,6)=dot(f(:,6),C(i,);
end
lny1=zeros(N,6);
for j=1:n
for i=1:N
if E(i,j)>0
lny1(i,j)=A(i,7)*(1-log(E(i,j))-F(i,j));
else
lny1(i,j)=0;
end
end
end
for i=1:N
lny1(i,6)=A(i,7)*(1-log(E(i,6))-F(i,6));
end
lny=zeros(1,n);
for j=1:n
lny(j)=lny2(j)+dot(A(1:N,j),lny1(:,6))-dot(A(1:N,j),lny1(:,j));
end
yc=zeros(1,n);
for i=1:n
yc(i)=exp(lny(i));
end
f=yc;
%---------------------------------------------------------------------------------------
function f=Temper(xx,t,n,Pt,R,Antoine,A,B,D,N)
t1=dot(xx,t);
yc=ActivityCoefficient(xx,t1,n,R,A,B,D,N);
Q=ones(1,n);
VPsat=VaporPressure(t1,n,Antoine);
VPsat(1)=Pt*VPsat(1)/dot(dot(dot(yc,xx),VPsat),1./Q);
T=Antoine(2,1)/(Antoine(1,1)-lg(100000*VPsat(1)))-Antoine(3,1);
f=T;