问题描述
请帮助我解决这个问题，因为我需要使用 nlinfit 函数来拟合未知参数，但它显示了一些错误。虽然昨天我得到了一些要拟合的参数值，但现在如果我正在运行它，函数输出将出现一些问题，用于拟合 NaN 答案仅用于最后一次迭代。 X数据是一个三列串联矩阵作为自变量，yk是因变量，taua是待拟合参数个数的初始估计矩阵。
function [yk]=activity_coefficientE(taua,x)
T=523;
alpha12=0.3; alpha13=0.3; alpha21=0.3; alpha23=0.3; alpha31=0.3; alpha32=0.3;
alpha18=0.2; alpha81=0.2; alpha28=0.2; alpha82=0.2; alpha38=0.2; alpha83=0.3;
alpha19=0.2; alpha91=0.2; alpha29=0.2; alpha92=0.2; alpha39=0.2; alpha93=0.2;
alpha110=0.2;alpha101=0.2;alpha210=0.2;alpha102=0.2;alpha310=0.2;alpha103=0.2;
alpha113=0.2;alpha131=0.2;alpha213=0.2;alpha132=0.2;alpha313=0.2;alpha133=0.2;
alpha114=0.2;alpha141=0.2;alpha214=0.2;alpha142=0.2;alpha314=0.2;alpha143=0.2;
alpha115=0.2;alpha151=0.2;alpha215=0.2;alpha152=0.2;alpha315=0.2;alpha153=0.2;
alpha117=0.2;alpha171=0.2;alpha217=0.2;alpha172=0.2;alpha317=0.2;alpha173=0.2;
alpha118=0.2;alpha181=0.2;alpha218=0.2;alpha182=0.2;alpha318=0.2;alpha183=0.2;
alpha810=0.2;alpha915=0.2;alpha1314=0.2;alpha108=0.2;alpha159=0.2;alpha1413=0.2;
alpha1718=0.2;alpha1817=0.2;
tau12=0; tau13=0; tau21=0; tau23=0; tau31=0; tau32=0;
%taua=randi([-5,5],1,112)
tau18=taua(1)+taua(57)/T;
tau81=taua(2)+taua(58)/T;
tau28=taua(3)+taua(59)/T;
tau82=taua(4)+taua(60)/T;
tau38=taua(5)+taua(61)/T;
tau83=taua(6)+taua(62)/T;
tau19=taua(7)+taua(63)/T;
tau91=taua(8)+taua(64)/T;
tau29=taua(9)+taua(65)/T;
tau92=taua(10)+taua(66)/T;
tau39=taua(11)+taua(67)/T;
tau93=taua(12)+taua(68)/T;
tau110=taua(13)+taua(69)/T;
tau101=taua(14)+taua(70)/T;
tau210=taua(15)+taua(71)/T;
tau102=taua(16)+taua(72)/T;
tau310=taua(17)+taua(73)/T;
tau103=taua(18)+taua(74)/T;
tau113=taua(19)+taua(75)/T;
tau131=taua(20)+taua(76)/T;
tau213=taua(21)+taua(77)/T;
tau132=taua(22)+taua(78)/T;
tau313=taua(23)+taua(79)/T;
tau133=taua(24)+taua(80)/T;
tau114=taua(25)+taua(81)/T;
tau141=taua(26)+taua(82)/T;
tau214=taua(27)+taua(83)/T;
tau142=taua(28)+taua(84)/T;
tau314=taua(29)+taua(85)/T;
tau143=taua(30)+taua(86)/T;
tau115=taua(31)+taua(87)/T;
tau151=taua(32)+taua(88)/T;
tau215=taua(33)+taua(89)/T;
tau152=taua(34)+taua(90)/T;
tau315=taua(35)+taua(91)/T;
tau153=taua(36)+taua(92)/T;
tau117=taua(37)+taua(93)/T;
tau171=taua(38)+taua(94)/T;
tau217=taua(39)+taua(95)/T;
tau172=taua(40)+taua(96)/T;
tau317=taua(41)+taua(97)/T;
tau173=taua(42)+taua(98)/T;
tau118=taua(43)+taua(99)/T;
tau181=taua(44)+taua(100)/T;
tau218=taua(45)+taua(101)/T;
tau182=taua(46)+taua(102)/T;
tau318=taua(47)+taua(103)/T;
tau183=taua(48)+taua(104)/T;
tau810=taua(49)+taua(105)/T;
tau108=taua(50)+taua(106)/T;
tau915=taua(51)+taua(107)/T;
tau159=taua(52)+taua(108)/T;
tau1314=taua(53)+taua(109)/T;
tau1413=taua(54)+taua(110)/T;
tau1718=taua(55)+taua(111)/T;
tau1817=taua(56)+taua(112)/T;
G12=exp(-(tau12*alpha12));
G21=exp(-(tau21*alpha21));
G13=exp(-(tau13*alpha13));
G31=exp(-(tau31*alpha31));
G23=exp(-(tau23*alpha23));
G32=exp(-(tau32*alpha32));
G18=exp(-(tau18*alpha18));
G81=exp(-(tau81*alpha81));
G28=exp(-(tau28*alpha28));
G82=exp(-(tau82*alpha82));
G38=exp(-(tau38*alpha83));
G83=exp(-(tau83*alpha83));
G19=exp(-(tau19*alpha19));
G91=exp(-(tau91*alpha91));
G29=exp(-(tau29*alpha29));
G92=exp(-(tau92*alpha92));
G39=exp(-(tau39*alpha39));
G93=exp(-(tau93*alpha93));
G110=exp(-(tau110*alpha110));
G101=exp(-(tau101*alpha101));
G210=exp(-(tau210*alpha210));
G102=exp(-(tau102*alpha102));
G310=exp(-(tau310*alpha310));
G103=exp(-(tau103*alpha103));
G113=exp(-(tau113*alpha113));
G131=exp(-(tau131*alpha131));
G213=exp(-(tau213*alpha213));
G132=exp(-(tau132*alpha132));
G313=exp(-(tau313*alpha313));
G133=exp(-(tau133*alpha133));
G114=exp(-(tau114*alpha114));
G141=exp(-(tau141*alpha141));
G214=exp(-(tau214*alpha214));
G142=exp(-(tau142*alpha142));
G314=exp(-(tau314*alpha314));
G143=exp(-(tau143*alpha143));
G115=exp(-(tau115*alpha115));
G151=exp(-(tau151*alpha151));
G215=exp(-(tau215*alpha215));
G152=exp(-(tau152*alpha152));
G315=exp(-(tau315*alpha315));
G153=exp(-(tau153*alpha153));
G117=exp(-(tau117*alpha117));
G171=exp(-(tau171*alpha171));
G217=exp(-(tau217*alpha217));
G172=exp(-(tau172*alpha172));
G317=exp(-(tau317*alpha317));
G173=exp(-(tau173*alpha173));
G118=exp(-(tau118*alpha118));
G181=exp(-(tau181*alpha181));
G218=exp(-(tau218*alpha218));
G182=exp(-(tau182*alpha182));
G318=exp(-(tau318*alpha318));
G183=exp(-(tau183*alpha183));
G810=exp(-(tau810*alpha810));
G108=exp(-(tau108*alpha108));
G915=exp(-(tau915*alpha915));
G159=exp(-(tau159*alpha159));
G1314=exp(-(tau1314*alpha1314));
G1413=exp(-(tau1413*alpha1413));
G1718=exp(-(tau1718*alpha1718));
G1817=exp(-(tau1817*alpha1817));
%calculating mole fractions of ionic species
x1=x(:,1);
x2=x(:,2);
x3=x(:,3);
%x1=[0.1577 0.1492 0.1462 0.1366 0.1299 0.1180 0.0863 0.0761 0.0550 ];
%x2=[0.8278 0.7945 0.7678 0.7450 0.6979 0.6309 0.4611 0.4114 0.2952 ];
%x3=[0.0145 0.0563 0.0860 0.1184 0.1722 0.2511 0.4526 0.5125 0.6498 ];
A=[0.0674243 0.0773881 0.0843400 0.0865343 0.0899223 0.0882858 0.0715087 0.0643867 0.0483658];
B=[0.0141081 0.0479814 0.0643151 0.0737477 0.0820756 0.0838701 0.0701576 0.0634457 0.0479639];
C=[0.0565665 0.0450072 0.0387724 0.0313828 0.02506094 0.0186280 0.0092734 0.0073438 0.0041595 ];
D=[0.0336447 0.0267694 0.0230611 0.0186659 0.0149058 0.0110795 0.0055157 0.0043679 0.0024739 ];
E=[0.0008148 0.0008756 0.00087131 0.0008794 0.0008711 0.0008441 0.0007384 0.0006997 0.0005980 ];
N=length(A);
x1n=zeros(N,1);x2n=zeros(N,1);x3n=zeros(N,1);
X1=zeros(N,1);X2=zeros(N,1);X3=zeros(N,1);X4=zeros(N,1);X5=zeros(N,1);X6=zeros(N,1);X7=zeros(N,1);
X12=zeros(N,1);X16=zeros(N,1);
for i=1:N
x1n(i)=(x1(i)-A(i)-D(i)-2*E(i)-C(i)+3*B(i))
x2n(i)=(x2(i)-A(i)-C(i)-D(i))
x3n(i)=(x3(i)-B(i))
X1(i)=(x1n(i)/(x1n(i)+x2n(i)+x3n(i)+2*A(i)+4*B(i)+2*C(i)+2*D(i)+2*E(i)))
X2(i)=(x2n(i)/(x1n(i)+x2n(i)+x3n(i)+2*A(i)+4*B(i)+2*C(i)+2*D(i)+2*E(i)))
X3(i)=(x3n(i)/(x1n(i)+x2n(i)+x3n(i)+2*A(i)+4*B(i)+2*C(i)+2*D(i)+2*E(i)))
X4(i)=(A(i)+D(i)+E(i)/(x1n(i)+x2n(i)+x3n(i)+2*A(i)+4*B(i)+2*C(i)+2*D(i)+2*E(i)))
X5(i)=(C(i)/(x1n(i)+x2n(i)+x3n(i)+2*A(i)+4*B(i)+2*C(i)+2*D(i)+2*E(i)))
X6(i)=(A(i)-B(i)/(x1n(i)+x2n(i)+x3n(i)+2*A(i)+4*B(i)+2*C(i)+2*D(i)+2*E(i)))
X7(i)=(B(i)/(x1n(i)+x2n(i)+x3n(i)+2*A(i)+4*B(i)+2*C(i)+2*D(i)+2*E(i)))
X12(i)=(E(i)/(x1n(i)+x2n(i)+x3n(i)+2*A(i)+4*B(i)+2*C(i)+2*D(i)+2*E(i)))
X16(i)=(C(i)+D(i)/(x1n(i)+x2n(i)+x3n(i)+2*A(i)+4*B(i)+2*C(i)+2*D(i)+2*E(i)))
end
yc4=X4./(X4+X5);
yc5=X5./(X4+X5);
yc6=X6./(X6+X7+X12+X16);
yc7=X7./(X6+X7+X12+X16);
yc12=X12./(X6+X7+X12+X16);
yc16=X16./(X6+X7+X12+X16);
alpha14=yc6.*alpha18+yc7.*alpha19+yc12.*alpha113+yc16.*alpha117;
%alpha41=alpha14;
alpha24=yc6.*alpha28+yc7.*alpha29+yc12.*alpha213+yc16.*alpha217;
%alpha42=alpha24;
alpha34=yc6.*alpha38+yc7.*alpha39+yc12.*alpha313+yc16.*alpha317;
%alpha43=alpha34;
alpha15=yc6.*alpha110+yc7.*alpha115+yc12.*alpha114+yc16.*alpha118;
%alpha51=alpha15;
alpha25=yc6.*alpha210+yc7.*alpha215+yc12.*alpha214+yc16.*alpha218;
%alpha52=alpha25;
alpha35=yc6.*alpha310+yc7.*alpha315+yc12.*alpha314+yc16.*alpha318;
%alpha53=alpha35;
alpha16=yc4.*alpha81+yc5.*alpha101;
%alpha61=alpha16;
alpha26=yc4.*alpha82+yc5.*alpha102;
%alpha62=alpha26;
alpha36=yc4.*alpha83+yc5.*alpha103;
%alpha63=alpha36;
alpha17=yc4.*alpha91+yc5.*alpha151;
%alpha71=alpha17;
alpha27=yc4.*alpha92+yc5.*alpha152;
%alpha72=alpha27;
alpha37=yc4.*alpha93+yc5.*alpha153;
%alpha73=alpha37;
alpha112=yc4.*alpha131+yc5.*alpha141;
%alpha121=alpha112;
alpha212=yc4.*alpha132+yc5.*alpha142;
%alpha122=alpha212;
alpha312=yc4.*alpha133+yc5.*alpha143;
%alpha123=alpha312;
alpha116=yc4.*alpha171+yc5.*alpha181;
%alpha161=alpha116;
alpha216=yc4.*alpha172+yc5.*alpha182;
%alpha162=alpha216;
alpha316=yc4.*alpha173+yc5.*alpha183;
%alpha163=alpha316;
alpha46=yc5.*alpha810;
%alpha64=alpha46;
alpha47=yc5.*alpha915;
%alpha74=alpha47;
alpha412=yc5.*alpha1314;
%alpha124=alpha412;
alpha416=yc5.*alpha1718;
%alpha164=alpha416;
alpha56=yc4.*alpha108;
%alpha65=alpha56;
alpha57=yc4.*alpha159;
%alpha75=alpha57;
alpha512=yc4.*alpha1413;
%alpha125=alpha512;
alpha516=yc4.*alpha1817;
%alpha165=alpha516;
G14=yc6.*G18+yc7.*G19+yc12.*G113+yc16.*G117;
%G41=G14;
G24=yc6.*G28+yc7.*G29+yc12.*G213+yc16.*G217;
%G42=G24;
G34=yc6.*G38+yc7.*G39+yc12.*G313+yc16.*G317;
%G43=G34;
G15=yc6.*G110+yc7.*G115+yc12.*G114+yc16.*G118;
%G51=G15;
G25=yc6.*G210+yc7.*G215+yc12.*G214+yc16.*G218;
%G52=G25;
G35=yc6.*G310+yc7.*G315+yc12.*G314+yc16.*G318;
%G53=G35;
G16=yc4.*G81+yc5.*G101;
%G61=G16;
G26=yc4.*G82+yc5.*G102;
%G62=G26;
G36=yc4.*G83+yc5.*G103;
%G63=G36;
G17=yc4.*G91+yc5.*G151;
%G71=G17;
G27=yc4.*G92+yc5.*G152;
%G72=G27;
G37=yc4.*G93+yc5.*G153;
%G73=G37;
G112=yc4.*G131+yc5.*G141;
%G121=G112;
G212=yc4.*G132+yc5.*G142;
%G122=G212;
G312=yc4.*G133+yc5.*G143;
%G123=G312;
G116=yc4.*G171+yc5.*G181;
%G161=G116;
G216=yc4.*G172+yc5.*G182;
%G162=G216;
G316=yc4.*G173+yc5.*G183;
%G163=G316;
G46=yc5.*G810;
%G64=G46;
G47=yc5.*G915;
%G74=G47;
G412=yc5.*G1314;
%G124=G412;
G416=yc5.*G1718;
%G164=G416;
G56=yc4.*G108;
%G65=G56;
G57=yc4.*G159;
%G75=G57;
G512=yc4.*G1413;
%G125=G512;
G516=yc4.*G1817;
%G165=G516;
tau14=-log(G14)./alpha14;
%tau41=tau14;
tau24=-log(G24)./alpha24;
%tau42=tau24;
tau34=-log(G34)./alpha34;
%tau43=tau34;
tau15=-log(G15)./alpha15;
%tau51=tau15;
tau25=-log(G25)./alpha25;
%tau52=tau25;
tau35=-log(G35)./alpha35;
%tau53=tau35;
tau16=-log(G16)./alpha16;
%tau61=tau16;
tau26=-log(G26)./alpha26;
%tau62=tau26;
tau36=-log(G36)./alpha36;
%tau63=tau36;
tau17=-log(G17)./alpha17;
%tau71=tau17;
tau27=-log(G27)./alpha27;
%tau72=tau27;
tau37=-log(G37)./alpha37;
%tau73=tau37;
tau112=-log(G112)./alpha112;
%tau121=tau112;
tau212=-log(G212)./alpha212;
%tau122=tau212;
tau312=-log(G312)./alpha312;
%tau123=tau312;
tau116=-log(G116)./alpha116;
%tau161=tau116;
tau216=-log(G216)./alpha216;
%tau162=tau216;
tau316=-log(G316)./alpha316;
%tau163=tau316;
tau46=-log(G46)./alpha46;
%tau64=tau46;
tau47=-log(G47)./alpha47;
%tau74=tau47;
tau412=-log(G412)./alpha412;
%tau124=tau412;
tau416=-log(G416)./alpha416;
%tau164=tau416;
tau56=-log(G56)./alpha56;
%tau65=tau56;
tau57=-log(G57)./alpha57;
%tau75=tau57;
tau512=-log(G512)./alpha512;
%tau125=tau512;
tau516=-log(G516)./alpha516;
%tau165=tau516;
ln_y1_1=G12.*X2.*tau12+ G31.*X3.*tau13+ G14.*X4.*tau14+G15.*X5.*tau15+G16.*X6.*tau16+G17.*X7.*tau17+G112.*X12.*tau112+G116.*X16.*tau116;
ln_y1_2=G12.*X2+ G13.*X3+ G14.*X4+G15.*X5+G16.*X6+G17.*X7+G112.*X12+G116.*X16;
ln_y2_1=G21.*X1.*tau12+ G32.*X3.*tau32+ G24.*X4.*tau24+G25.*X5.*tau25+G26.*X6.*tau26+G27.*X7.*tau27+G212.*X12.*tau212+G216.*X16.*tau216;
ln_y2_2=G12.*X1+ G23.*X3+G24.*X4+G25.*X5+G26.*X6+G27.*X7+G212.*X12+G216.*X16;
ln_y3_1=G13.*X1.*tau13+ G23.*X3.*tau23+ G34.*X4.*tau34+G35.*X5.*tau35+G36.*X6.*tau36+G37.*X7.*tau37+G312.*X12.*tau312+G316.*X16.*tau316;
ln_y3_2=G13.*X1+ G23.*X3+ G34.*X4+G35.*X5+G36.*X6+G37.*X7+G312.*X12+G316.*X16;
ln_y4_1=G14.*X1.*tau14+G24.*X2.*tau24+G34.*X3.*tau34+G46.*X6.*tau46+G47.*X7.*tau47+G412.*X12.*tau412+G416.*X16.*tau416;
ln_y4_2=G14.*X1+G24.*X2+G34.*X3+G46.*X6+G47.*X7+G412.*X12+G416.*X16;
ln_y5_1=G15.*X1.*tau15+G25.*X2.*tau25+G35.*X3.*tau35+G56.*X6.*tau56+G57.*X7.*tau57+G512.*X12.*tau512+G516.*X16.*tau516;
ln_y5_2=G15.*X1+G25.*X2+G35.*X3+G56.*X6+G57.*X7+G512.*X12+G516.*X16;
ln_y6_1=G16.*X1.*tau16+G26.*X2.*tau26+G36.*X3.*tau36+G46.*X4.*tau46+G56.*X5.*tau56;
ln_y6_2=G16.*X1+G26.*X2+G36.*X3+G46.*X4+G56.*X5;
ln_y7_1=G17.*X1.*tau17+G27.*X2.*tau27+G37.*X3.*tau37+G47.*X4.*tau47+G57.*X5.*tau57;
ln_y7_2=G17.*X1+G27.*X2+G37.*X3+G47.*X4+G57.*X5;
ln_y12_1=G112.*X1.*tau112+G212.*X2.*tau212+G312.*X3.*tau312+G412.*X4.*tau412+G512.*X5.*tau512;
ln_y12_2=G112.*X1+G212.*X2+G312.*X3+G412.*X4+G512.*X5;
ln_y16_1=G116.*X1.*tau116+G216.*X2.*tau216+G316.*X3.*tau316+G416.*X4.*tau416+G516.*X5.*tau516;
ln_y16_2=G116.*X1+G216.*X2+G316.*X3+G416.*X4+G516.*X5;
ln_y1_3=(((X2.*G12)./ln_y2_2).*(tau12-(ln_y2_1)./(ln_y2_2)))+(((X3.*G13)./ln_y3_2).*(tau13-(ln_y3_1)./(ln_y3_2)));
ln_y1_4=(((X6.*G16)./ln_y6_2).*(tau16- (ln_y6_1./ln_y6_2))) + (((X7.*G17)./ln_y7_2).*(tau17- (ln_y7_1./ln_y7_2)))+(((X12.*G12)./ln_y12_2).*(tau112- (ln_y12_1./ln_y12_2)))+(((X16.*G16)./ln_y16_2).*(tau116- (ln_y16_1./ln_y16_2)));
ln_y1_5=(((X4.*G14)./ln_y4_2).*(tau14- (ln_y4_1./ln_y4_2))) + (((X5.*G15)./ln_y5_2).*(tau15- (ln_y5_1./ln_y5_2)));
yk=exp((ln_y1_1./ln_y1_2) + ln_y1_3 + ln_y1_4+ ln_y1_5) % activity coefficient for H2O
end
........................................
Another function where above function to be called.....
% calling the function act_coeff to estimate the binary interaction parameters
for i=1:112
filename = 'EagelsDATA.xlsx'; %reading VLE data from excel file
Data = xlsread(filename);
x(:,1) = Data([10:15 17:19],16);
x(:,2) = Data([10:15 17:19],1);
x(:,3)= Data([10:15 17:19],2);
taua=(randi([-5,112));
yk=[0.0606 (values calculated from above function and will be used for fitting)
0.4327
0.6545
0.9417
1.2570
1.6881
1.9108
1.7777
1.3821]
% taua =[ -2 3 4 -3 -4 1 4 -2 4 -4 -1 4 5 -3 3 2 -5 3 -4
% 1 4 1 5 -1 -1 -3 2 -3 4 3 4 2 5 4 -2 4 3 -1
% 1 0 -5 -5 -5 -3 4 2 1 4 0 2 -3 -4 5 0 -3 2 5
% 1 0 5 1 -3 5 4 1 5 2 3 2 0 -5 -4 -2 1 -2 5
%-5 5 -2 -2 4 1 -1 3 -1 1 5 -1 0 -1 4 5 5 1 4
% 1 0 4 -4 4 0 -1 -2 -5 -3 -4 -5
% -5 0 -2 0 -5] (random values for which yk was calculted from the command
taua= randi([-5,112))
try % try-catch used to continue the loop without stopping on encountering an error
[taua1]= nlinfit(x,yk,@activity_coefficientE,taua)
catch exception
continue
end
end
I am not able to attach excel sheet here so data from excel sheet is as:
x =[0.1577 0.1492 0.1462 0.1366 0.1299 0.1180 0.0863 0.0761 0.0550; column 1
0.8278 0.7945 0.7678 0.7450 0.6979 0.6309 0.4611 0.4114 0.2952 ; column 2
0.0145 0.0563 0.0860 0.1184 0.1722 0.2511 0.4526 0.5125 0.6498 ]; column 3
解决方法
我发现你所做的有 3 个主要问题。
问题 #1 - 错误 您收到错误的原因是您的函数“activity_coefficientE”有时会返回 NaN 或 inf 值。我的建议是寻找这些值并将“yk”的值设置为一个较大的值，以便“nlinfit”中的优化器将远离产生无限或 NaN 值的系数。这是函数底部的代码，可以避免崩溃：
if any(~isfinite(yk))
yk = 10 * ones(size(yk));
end
问题 2 - 随机初始猜测 使用随机数作为初始条件的麻烦在于，每次运行它时都会得到不同的答案，因此有时它有效，有时则无效。如果设置随机数生成器种子，则每次运行脚本时都可以获得相同的随机数。如果您更改种子，则可以获得一组不同的随机数。我将您的主要脚本缩短为这样，我尝试 100 个不同的随机种子（并存储每次尝试的结果）以查看结果：
for i=1:100
rng(i)
taua = randi([-5,5],1,112);
taua1(i,:) = nlinfit(x,yk,@activity_coefficientE,taua);
end
“taua1”的每一行是一组 111 个系数。
问题 #3 - 尝试用 112 个系数拟合 9 个点 每次调用 nlinfit 时，您都会收到此警告：
警告：排名不足
因为您要求 nlinfit 找到的系数 (112) 多于您拟合的数据点 (9)。这就像试图找到最适合 2 个点的二阶方程，有无数个解。曲线拟合时，您应该拥有比系数更多的数据点，以确保您没有拟合噪声。您需要“yk”和“x”中的更多数据点和/或更少的系数来拟合。我做了很多曲线拟合，但从未见过有 112 个系数的方程，所以我认为您没有正确解决问题。也许这 112 个系数并不是真正独立的，或者您想要找到 112 个数据点和 9 个系数。
为了完整起见，这是我为处理此解决方案而创建的 activity_coefficientE.m 函数的编辑版本。一般来说，我从未见过有这么多名称相似的变量的 Matlab 代码。通过使用向量运算可以大大简化大部分代码。我的大部分更改都涉及格式化、添加分号以及最后检查非有限值。
function yk=activity_coefficientE(taua,x)
T=523;
alpha12=0.3; alpha13=0.3; alpha21=0.3; alpha23=0.3; alpha31=0.3; alpha32=0.3;
alpha18=0.2; alpha81=0.2; alpha28=0.2; alpha82=0.2; alpha38=0.2; alpha83=0.3;
alpha19=0.2; alpha91=0.2; alpha29=0.2; alpha92=0.2; alpha39=0.2; alpha93=0.2;
alpha110=0.2;alpha101=0.2;alpha210=0.2;alpha102=0.2;alpha310=0.2;alpha103=0.2;
alpha113=0.2;alpha131=0.2;alpha213=0.2;alpha132=0.2;alpha313=0.2;alpha133=0.2;
alpha114=0.2;alpha141=0.2;alpha214=0.2;alpha142=0.2;alpha314=0.2;alpha143=0.2;
alpha115=0.2;alpha151=0.2;alpha215=0.2;alpha152=0.2;alpha315=0.2;alpha153=0.2;
alpha117=0.2;alpha171=0.2;alpha217=0.2;alpha172=0.2;alpha317=0.2;alpha173=0.2;
alpha118=0.2;alpha181=0.2;alpha218=0.2;alpha182=0.2;alpha318=0.2;alpha183=0.2;
alpha810=0.2;alpha915=0.2;alpha1314=0.2;alpha108=0.2;alpha159=0.2;alpha1413=0.2;
alpha1718=0.2;alpha1817=0.2;
tau12=0; tau13=0; tau21=0; tau23=0; tau31=0; tau32=0;
tau18=taua(1)+taua(57)/T;
tau81=taua(2)+taua(58)/T;
tau28=taua(3)+taua(59)/T;
tau82=taua(4)+taua(60)/T;
tau38=taua(5)+taua(61)/T;
tau83=taua(6)+taua(62)/T;
tau19=taua(7)+taua(63)/T;
tau91=taua(8)+taua(64)/T;
tau29=taua(9)+taua(65)/T;
tau92=taua(10)+taua(66)/T;
tau39=taua(11)+taua(67)/T;
tau93=taua(12)+taua(68)/T;
tau110=taua(13)+taua(69)/T;
tau101=taua(14)+taua(70)/T;
tau210=taua(15)+taua(71)/T;
tau102=taua(16)+taua(72)/T;
tau310=taua(17)+taua(73)/T;
tau103=taua(18)+taua(74)/T;
tau113=taua(19)+taua(75)/T;
tau131=taua(20)+taua(76)/T;
tau213=taua(21)+taua(77)/T;
tau132=taua(22)+taua(78)/T;
tau313=taua(23)+taua(79)/T;
tau133=taua(24)+taua(80)/T;
tau114=taua(25)+taua(81)/T;
tau141=taua(26)+taua(82)/T;
tau214=taua(27)+taua(83)/T;
tau142=taua(28)+taua(84)/T;
tau314=taua(29)+taua(85)/T;
tau143=taua(30)+taua(86)/T;
tau115=taua(31)+taua(87)/T;
tau151=taua(32)+taua(88)/T;
tau215=taua(33)+taua(89)/T;
tau152=taua(34)+taua(90)/T;
tau315=taua(35)+taua(91)/T;
tau153=taua(36)+taua(92)/T;
tau117=taua(37)+taua(93)/T;
tau171=taua(38)+taua(94)/T;
tau217=taua(39)+taua(95)/T;
tau172=taua(40)+taua(96)/T;
tau317=taua(41)+taua(97)/T;
tau173=taua(42)+taua(98)/T;
tau118=taua(43)+taua(99)/T;
tau181=taua(44)+taua(100)/T;
tau218=taua(45)+taua(101)/T;
tau182=taua(46)+taua(102)/T;
tau318=taua(47)+taua(103)/T;
tau183=taua(48)+taua(104)/T;
tau810=taua(49)+taua(105)/T;
tau108=taua(50)+taua(106)/T;
tau915=taua(51)+taua(107)/T;
tau159=taua(52)+taua(108)/T;
tau1314=taua(53)+taua(109)/T;
tau1413=taua(54)+taua(110)/T;
tau1718=taua(55)+taua(111)/T;
tau1817=taua(56)+taua(112)/T;
G12=exp(-(tau12*alpha12));
G21=exp(-(tau21*alpha21));
G13=exp(-(tau13*alpha13));
G31=exp(-(tau31*alpha31));
G23=exp(-(tau23*alpha23));
G32=exp(-(tau32*alpha32));
G18=exp(-(tau18*alpha18));
G81=exp(-(tau81*alpha81));
G28=exp(-(tau28*alpha28));
G82=exp(-(tau82*alpha82));
G38=exp(-(tau38*alpha83));
G83=exp(-(tau83*alpha83));
G19=exp(-(tau19*alpha19));
G91=exp(-(tau91*alpha91));
G29=exp(-(tau29*alpha29));
G92=exp(-(tau92*alpha92));
G39=exp(-(tau39*alpha39));
G93=exp(-(tau93*alpha93));
G110=exp(-(tau110*alpha110));
G101=exp(-(tau101*alpha101));
G210=exp(-(tau210*alpha210));
G102=exp(-(tau102*alpha102));
G310=exp(-(tau310*alpha310));
G103=exp(-(tau103*alpha103));
G113=exp(-(tau113*alpha113));
G131=exp(-(tau131*alpha131));
G213=exp(-(tau213*alpha213));
G132=exp(-(tau132*alpha132));
G313=exp(-(tau313*alpha313));
G133=exp(-(tau133*alpha133));
G114=exp(-(tau114*alpha114));
G141=exp(-(tau141*alpha141));
G214=exp(-(tau214*alpha214));
G142=exp(-(tau142*alpha142));
G314=exp(-(tau314*alpha314));
G143=exp(-(tau143*alpha143));
G115=exp(-(tau115*alpha115));
G151=exp(-(tau151*alpha151));
G215=exp(-(tau215*alpha215));
G152=exp(-(tau152*alpha152));
G315=exp(-(tau315*alpha315));
G153=exp(-(tau153*alpha153));
G117=exp(-(tau117*alpha117));
G171=exp(-(tau171*alpha171));
G217=exp(-(tau217*alpha217));
G172=exp(-(tau172*alpha172));
G317=exp(-(tau317*alpha317));
G173=exp(-(tau173*alpha173));
G118=exp(-(tau118*alpha118));
G181=exp(-(tau181*alpha181));
G218=exp(-(tau218*alpha218));
G182=exp(-(tau182*alpha182));
G318=exp(-(tau318*alpha318));
G183=exp(-(tau183*alpha183));
G810=exp(-(tau810*alpha810));
G108=exp(-(tau108*alpha108));
G915=exp(-(tau915*alpha915));
G159=exp(-(tau159*alpha159));
G1314=exp(-(tau1314*alpha1314));
G1413=exp(-(tau1413*alpha1413));
G1718=exp(-(tau1718*alpha1718));
G1817=exp(-(tau1817*alpha1817));
%calculating mole fractions of ionic species
x1=x(:,1);
x2=x(:,2);
x3=x(:,3);
A=[0.0674243 0.0773881 0.0843400 0.0865343 0.0899223 0.0882858 0.0715087 0.0643867 0.0483658];
B=[0.0141081 0.0479814 0.0643151 0.0737477 0.0820756 0.0838701 0.0701576 0.0634457 0.0479639];
C=[0.0565665 0.0450072 0.0387724 0.0313828 0.02506094 0.0186280 0.0092734 0.0073438 0.0041595 ];
D=[0.0336447 0.0267694 0.0230611 0.0186659 0.0149058 0.0110795 0.0055157 0.0043679 0.0024739 ];
E=[0.0008148 0.0008756 0.00087131 0.0008794 0.0008711 0.0008441 0.0007384 0.0006997 0.0005980 ];
N=length(A);
x1n=zeros(N,1);x2n=zeros(N,1);x3n=zeros(N,1);
X1=zeros(N,1);X2=zeros(N,1);X3=zeros(N,1);X4=zeros(N,1);X5=zeros(N,1);X6=zeros(N,1);X7=zeros(N,1);
X12=zeros(N,1);X16=zeros(N,1);
for i=1:N
x1n(i)=(x1(i)-A(i)-D(i)-2*E(i)-C(i)+3*B(i));
x2n(i)=(x2(i)-A(i)-C(i)-D(i));
x3n(i)=(x3(i)-B(i));
X1(i)=(x1n(i)/(x1n(i)+x2n(i)+x3n(i)+2*A(i)+4*B(i)+2*C(i)+2*D(i)+2*E(i)));
X2(i)=(x2n(i)/(x1n(i)+x2n(i)+x3n(i)+2*A(i)+4*B(i)+2*C(i)+2*D(i)+2*E(i)));
X3(i)=(x3n(i)/(x1n(i)+x2n(i)+x3n(i)+2*A(i)+4*B(i)+2*C(i)+2*D(i)+2*E(i)));
X4(i)=(A(i)+D(i)+E(i)/(x1n(i)+x2n(i)+x3n(i)+2*A(i)+4*B(i)+2*C(i)+2*D(i)+2*E(i)));
X5(i)=(C(i)/(x1n(i)+x2n(i)+x3n(i)+2*A(i)+4*B(i)+2*C(i)+2*D(i)+2*E(i)));
X6(i)=(A(i)-B(i)/(x1n(i)+x2n(i)+x3n(i)+2*A(i)+4*B(i)+2*C(i)+2*D(i)+2*E(i)));
X7(i)=(B(i)/(x1n(i)+x2n(i)+x3n(i)+2*A(i)+4*B(i)+2*C(i)+2*D(i)+2*E(i)));
X12(i)=(E(i)/(x1n(i)+x2n(i)+x3n(i)+2*A(i)+4*B(i)+2*C(i)+2*D(i)+2*E(i)));
X16(i)=(C(i)+D(i)/(x1n(i)+x2n(i)+x3n(i)+2*A(i)+4*B(i)+2*C(i)+2*D(i)+2*E(i)));
end
yc4=X4./(X4+X5);
yc5=X5./(X4+X5);
yc6=X6./(X6+X7+X12+X16);
yc7=X7./(X6+X7+X12+X16);
yc12=X12./(X6+X7+X12+X16);
yc16=X16./(X6+X7+X12+X16);
alpha14=yc6.*alpha18+yc7.*alpha19+yc12.*alpha113+yc16.*alpha117;
alpha24=yc6.*alpha28+yc7.*alpha29+yc12.*alpha213+yc16.*alpha217;
alpha34=yc6.*alpha38+yc7.*alpha39+yc12.*alpha313+yc16.*alpha317;
alpha15=yc6.*alpha110+yc7.*alpha115+yc12.*alpha114+yc16.*alpha118;
alpha25=yc6.*alpha210+yc7.*alpha215+yc12.*alpha214+yc16.*alpha218;
alpha35=yc6.*alpha310+yc7.*alpha315+yc12.*alpha314+yc16.*alpha318;
alpha16=yc4.*alpha81+yc5.*alpha101;
alpha26=yc4.*alpha82+yc5.*alpha102;
alpha36=yc4.*alpha83+yc5.*alpha103;
alpha17=yc4.*alpha91+yc5.*alpha151;
alpha27=yc4.*alpha92+yc5.*alpha152;
alpha37=yc4.*alpha93+yc5.*alpha153;
alpha112=yc4.*alpha131+yc5.*alpha141;
alpha212=yc4.*alpha132+yc5.*alpha142;
alpha312=yc4.*alpha133+yc5.*alpha143;
alpha116=yc4.*alpha171+yc5.*alpha181;
alpha216=yc4.*alpha172+yc5.*alpha182;
alpha316=yc4.*alpha173+yc5.*alpha183;
alpha46=yc5.*alpha810;
alpha47=yc5.*alpha915;
alpha412=yc5.*alpha1314;
alpha416=yc5.*alpha1718;
alpha56=yc4.*alpha108;
alpha57=yc4.*alpha159;
alpha512=yc4.*alpha1413;
alpha516=yc4.*alpha1817;
G14=yc6.*G18+yc7.*G19+yc12.*G113+yc16.*G117;
G24=yc6.*G28+yc7.*G29+yc12.*G213+yc16.*G217;
G34=yc6.*G38+yc7.*G39+yc12.*G313+yc16.*G317;
G15=yc6.*G110+yc7.*G115+yc12.*G114+yc16.*G118;
G25=yc6.*G210+yc7.*G215+yc12.*G214+yc16.*G218;
G35=yc6.*G310+yc7.*G315+yc12.*G314+yc16.*G318;
G16=yc4.*G81+yc5.*G101;
G26=yc4.*G82+yc5.*G102;
G36=yc4.*G83+yc5.*G103;
G17=yc4.*G91+yc5.*G151;
G27=yc4.*G92+yc5.*G152;
G37=yc4.*G93+yc5.*G153;
G112=yc4.*G131+yc5.*G141;
G212=yc4.*G132+yc5.*G142;
G312=yc4.*G133+yc5.*G143;
G116=yc4.*G171+yc5.*G181;
G216=yc4.*G172+yc5.*G182;
G316=yc4.*G173+yc5.*G183;
G46=yc5.*G810;
G47=yc5.*G915;
G412=yc5.*G1314;
G416=yc5.*G1718;
G56=yc4.*G108;
G57=yc4.*G159;
G512=yc4.*G1413;
G516=yc4.*G1817;
tau14=-log(G14)./alpha14;
tau24=-log(G24)./alpha24;
tau34=-log(G34)./alpha34;
tau15=-log(G15)./alpha15;
tau25=-log(G25)./alpha25;
tau35=-log(G35)./alpha35;
tau16=-log(G16)./alpha16;
tau26=-log(G26)./alpha26;
tau36=-log(G36)./alpha36;
tau17=-log(G17)./alpha17;
tau27=-log(G27)./alpha27;
tau37=-log(G37)./alpha37;
tau112=-log(G112)./alpha112;
tau212=-log(G212)./alpha212;
tau312=-log(G312)./alpha312;
tau116=-log(G116)./alpha116;
tau216=-log(G216)./alpha216;
tau316=-log(G316)./alpha316;
tau46=-log(G46)./alpha46;
tau47=-log(G47)./alpha47;
tau412=-log(G412)./alpha412;
tau416=-log(G416)./alpha416;
tau56=-log(G56)./alpha56;
tau57=-log(G57)./alpha57;
tau512=-log(G512)./alpha512;
tau516=-log(G516)./alpha516;
ln_y1_1=G12.*X2.*tau12+ G31.*X3.*tau13+ G14.*X4.*tau14+G15.*X5.*tau15+G16.*X6.*tau16+G17.*X7.*tau17+G112.*X12.*tau112+G116.*X16.*tau116;
ln_y1_2=G12.*X2+ G13.*X3+ G14.*X4+G15.*X5+G16.*X6+G17.*X7+G112.*X12+G116.*X16;
ln_y2_1=G21.*X1.*tau12+ G32.*X3.*tau32+ G24.*X4.*tau24+G25.*X5.*tau25+G26.*X6.*tau26+G27.*X7.*tau27+G212.*X12.*tau212+G216.*X16.*tau216;
ln_y2_2=G12.*X1+ G23.*X3+G24.*X4+G25.*X5+G26.*X6+G27.*X7+G212.*X12+G216.*X16;
ln_y3_1=G13.*X1.*tau13+ G23.*X3.*tau23+ G34.*X4.*tau34+G35.*X5.*tau35+G36.*X6.*tau36+G37.*X7.*tau37+G312.*X12.*tau312+G316.*X16.*tau316;
ln_y3_2=G13.*X1+ G23.*X3+ G34.*X4+G35.*X5+G36.*X6+G37.*X7+G312.*X12+G316.*X16;
ln_y4_1=G14.*X1.*tau14+G24.*X2.*tau24+G34.*X3.*tau34+G46.*X6.*tau46+G47.*X7.*tau47+G412.*X12.*tau412+G416.*X16.*tau416;
ln_y4_2=G14.*X1+G24.*X2+G34.*X3+G46.*X6+G47.*X7+G412.*X12+G416.*X16;
ln_y5_1=G15.*X1.*tau15+G25.*X2.*tau25+G35.*X3.*tau35+G56.*X6.*tau56+G57.*X7.*tau57+G512.*X12.*tau512+G516.*X16.*tau516;
ln_y5_2=G15.*X1+G25.*X2+G35.*X3+G56.*X6+G57.*X7+G512.*X12+G516.*X16;
ln_y6_1=G16.*X1.*tau16+G26.*X2.*tau26+G36.*X3.*tau36+G46.*X4.*tau46+G56.*X5.*tau56;
ln_y6_2=G16.*X1+G26.*X2+G36.*X3+G46.*X4+G56.*X5;
ln_y7_1=G17.*X1.*tau17+G27.*X2.*tau27+G37.*X3.*tau37+G47.*X4.*tau47+G57.*X5.*tau57;
ln_y7_2=G17.*X1+G27.*X2+G37.*X3+G47.*X4+G57.*X5;
ln_y12_1=G112.*X1.*tau112+G212.*X2.*tau212+G312.*X3.*tau312+G412.*X4.*tau412+G512.*X5.*tau512;
ln_y12_2=G112.*X1+G212.*X2+G312.*X3+G412.*X4+G512.*X5;
ln_y16_1=G116.*X1.*tau116+G216.*X2.*tau216+G316.*X3.*tau316+G416.*X4.*tau416+G516.*X5.*tau516;
ln_y16_2=G116.*X1+G216.*X2+G316.*X3+G416.*X4+G516.*X5;
ln_y1_3=(((X2.*G12)./ln_y2_2).*(tau12-(ln_y2_1)./(ln_y2_2)))+(((X3.*G13)./ln_y3_2).*(tau13-(ln_y3_1)./(ln_y3_2)));
ln_y1_4=(((X6.*G16)./ln_y6_2).*(tau16- (ln_y6_1./ln_y6_2))) + (((X7.*G17)./ln_y7_2).*(tau17- (ln_y7_1./ln_y7_2)))+(((X12.*G12)./ln_y12_2).*(tau112- (ln_y12_1./ln_y12_2)))+(((X16.*G16)./ln_y16_2).*(tau116- (ln_y16_1./ln_y16_2)));
ln_y1_5=(((X4.*G14)./ln_y4_2).*(tau14- (ln_y4_1./ln_y4_2))) + (((X5.*G15)./ln_y5_2).*(tau15- (ln_y5_1./ln_y5_2)));
yk=exp((ln_y1_1./ln_y1_2) + ln_y1_3 + ln_y1_4+ ln_y1_5)'; % activity coefficient for H2O
if any(~isfinite(yk))
yk = 10 * ones(size(yk));
end