斑竹及各位大侠帮小弟看一个程序
一个非线性方程组，有三百多个用fsolve解出的的效果不好，恳请斑竹及各位大侠推荐其他的求解方法，或修改的东西
源程序见附页
function CSTR_644
%对不同的循环圈采取不同的排出流量
%漩涡流量根据CFD模拟结果计算得到
%循环圈内的循环流量认为大小不等
%同时考虑氧对反应的影响，
z1=zeros(16,4,6);
z1(:,:,5)=ones(16,4)*2;
z1(:,:,6)=zeros(16,4);
x0=z1(:);
%options=optimset;%('Display','off','Jacobian','on');
%options.Tolx=1e-10;
%options.Tolfun=1e-20;
%options=optimset('MaxfunEvals',,'Tolx',,'Tolfun',,'Jacobian','off');
y=fsolve(@CSTR644,x0);
m=reshape(abs(y),16,4,6)
a1=sum(sum(m(:,:,1)));
a2=sum(sum(m(:,:,2)));
a3=sum(sum(m(:,:,3)));
a4=sum(sum(m(:,:,4)));
a5=sum(sum(m(:,:,5)));
a8=sum(sum(m(:,:,6)));
%反应收率
a6=a5/(a1)
%杂质含量
a7=a4/(a1)
function E=CSTR644(x)
c=reshape(x,16,4,6);
%设定气含率为5%均匀
Q1=25.3974;%第一主体循环流量
Q2=19.483;%第二主体循环流量
Q3=18.6079;%thE(i,j,h) third circlE(i,j,h) flow
Q4=18.9814;%thE(i,j,h) forth circlE(i,j,h) flow
F=0.005457;%thE(i,j,h) matE(i,j,h)rial fE(i,j,h)E(i,j,h)d and output (kg/s)
V1=4.684/64;%单元格体积
beta=0.2;%交互系数
c0=2.180;%加料中PX的初始浓度相对于醋酸mole/kg
b=V1*380*0.2/37.98;%定义为单元格内氧气的浓度(nE(i,j,h)E(i,j,h)d to modify)
V=V1*0.95;%单元格液相
M=V*1050;
T=463;% K
%k10=2.86e-3;%指前因子
k11=0.59+0.24*b-0.036*b.^2;%考虑氧浓度对反应的影响
%k(1)=k10*k11*exp(15.12-(17.82e+3)/T);%速率常数
k(1)=0.1716*k11/60;
%k20=1.16e-2;
k21=0.53+0.29*b-0.047*b.^2;
%k(2)=k20*k21*exp(13.20-(6.199e+3)/T);
k(2)=k21*0.7002/60;
%k30=5.88e-4;
k31=0.55+0.26*b-0.038*b.^2;
%k(3)=k30*k31*exp(18.83-(10.23e+3)/T);
k(3)=k31*0.0353/60;
%k40=5.49e-4;
k41=0.46+0.13*b-0.007*b.^2;
%k(4)=k40*k41*exp(20.41-(9.412e+3)/T);
k(4)=k41*0.3296/60;
%k50=6.44e-4;
k51=0.46+0.33*b-0.049*b.^2;
%k(5)=k50*k51*exp(15.10-(7.851e+3)/T);
k(5)=k51*0.0386/60;
%k60=2.14e-2;
k61=0.69+0.2*b-0.036*b.^2;
%k(6)=k60*k61*exp(15.51-(6.904e+3)/T);
k(6)=k61*1.2845/60;
%反应速率
j=linspace(1,4,4);
i=linspace(1,16,16);
r(i,j,1)=k(1)*c(i,j,1)./(1.4247*c(i,j,1)+4.8419*c(i,j,4)+0.0146);
%--------------------------------------------------------
r(i,j,2)=k(2)*c(i,j,2)./(1.4247*c(i,j,1)+4.8419*c(i,j,4)+0.0146).^0.5254;
%-----------------------------------------------------------
r(i,j,3)=k(3)*c(i,j,3)./(1.4247*c(i,j,1)+4.8419*c(i,j,4)+0.0146);
%-----------------------------------------------------
r(i,j,4)=k(4)*c(i,j,4)./(1.4247*c(i,j,1)+4.8419*c(i,j,4)+0.0146).^0.8111;
%----------------------------------------------------------
r(i,j,5)=k(5)*c(i,j,1)./(1.4247*c(i,j,1)+4.8419*c(i,j,4)+0.0146);
r(i,j,6)=k(6)*c(i,j,6)./(1.4247*c(i,j,1)+4.8419*c(i,j,4)+0.0146).^0.9302;
%----------------------------------------------------------------
%----------------------------------------------------------------
%balancE(i,j,h) E(i,j,h)quation of E(i,j,h)ach cE(i,j,h)ll
h=1;
j=1;
i=1;
E(i,j,h)=Q1*(c(i,j+1,h)-c(i,j,h))-M*r(i,j,1);
i=linspace(2,3,2);
E(i,j,h)=Q1*(c(i-1,j,h)-c(i,j,h))+beta*Q1*(c(i,j+1,h)-c(i,j,h))-M*r(i,j,1);
i=4;
E(i,j,h)=Q1*c(i-1,j,h)+F*c(i+1,j,h)-(Q1+F)*c(i,j,h)+10*beta*Q1*(c(i+1,j,h)-c(i,j,h))-M*r(i,j,1);
i=5;
E(i,j,h)=(Q2+2/3*F)*(c(i+1,j,h)-c(i,j,h))+10*beta*Q1*(c(i-1,j,h)-c(i,j,h))-M*r(i,j,1);
i=linspace(6,7,2);
E(i,j,h)=(Q2+2/3*F)*(c(i+1,j,h)-c(i,j,h))+beta*Q2*(c(i,j+1,h)-c(i,j,h))-M*r(i,j,1);
i=8;
E(i,j,h)=(Q2+2/3*F)*(c(i,j+1,h)-c(i,j,h))+beta*Q2*(c(i+1,j,h)-c(i,j,h))-M*r(i,j,1);
%第二个循环与第三个循环之间物料交换采用不同的流量以上一个循环为标准
i=9;
E(i,j,h)=Q3*(c(i,j+1,h)-c(i,j,h))+beta*Q2*(c(i-1,j,h)-c(i,j,h))-M*r(i,j,1);
i=linspace(10,11,2);
E(i,j,h)=Q3*(c(i-1,j,h)-c(i,j,h))+beta*Q3*(c(i,j+1,h)-c(i,j,h))-M*r(i,j,1);
i=12;
E(i,j,h)=Q3*c(i-1,j,h)+F*c(i+1,j,h)-(Q3+F)*c(i,j,h)+10*beta*Q3*(c(i+1,j,h)-c(i,j,h))-M*r(i,j,1);
i=13;
E(i,j,h)=(Q4+2/3*F)*(c(i+1,j,h)-c(i,j,h))+10*beta*Q3*(c(i-1,j,h)-c(i,j,h))-M*r(i,j,1);
i=linspace(14,15,2);
E(i,j,h)=(Q4+2/3*F)*(c(i+1,j,h)-c(i,j,h))+beta*Q4*(c(i,j+1,h)-c(i,j,h))-M*r(i,j,1);
i=16;
E(i,j,h)=(Q4+2/3*F)*(c(i,j+1,h)-c(i,j,h))-M*r(i,j,1);
%------------------------------------------------------------------------
j=2;
i=1;
E(i,j,h)=Q1*(c(i,j+1,h)-c(i,j,h))+beta*Q1*(c(i+1,j,h)-c(i,j,h))-M*r(i,j,1);
i=2;
E(i,j,h)=Q1*(c(i,j+1,h)-c(i,j,h))+beta*Q1*(c(i-1,j,h)+c(i,j-1,h)-2*c(i,j,h))-M*r(i,j,1);
i=3;
E(i,j,h)=Q1*c(i-1,j,h)+1/3*F*c(i+1,j,h)-(Q1+1/3*F)*c(i,j,h)+beta*Q1*(c(i,j-1,h)+c(i+1,j,h)-2*c(i,j,h))-M*r(i,j,1);
i=4;
E(i,j,h)=(Q1+F)*(c(i,j-1,h)-c(i,j,h))+10*beta*Q1*(c(i+1,j,h)-c(i,j,h))+beta*Q1*(c(i-1,j,h)-c(i,j,h))-M*r(i,j,1);
i=5;
E(i,j,h)=(Q2-1/3*F)*c(i,j-1,h)+1/3*F*c(i+1,j,h)-Q2*c(i,j,h)+10*beta*Q1*(c(i-1,j,h)-c(i,j,h))+beta*Q2*(c(i+1,j,h)-c(i,j,h))-M*r(i,j,1);
i=6;
E(i,j,h)=(Q2+1/3*F)*(c(i+1,j,h)-c(i,j,h))+beta*Q2*(c(i,j-1,h)+c(i-1,j,h)-2*c(i,j,h))-M*r(i,j,1);
i=7;
E(i,j,h)=(Q2+1/3*F)*(c(i,j+1,h)-c(i,j,h))+beta*Q2*(c(i,j-1,h)+c(i+1,j,h)-2*c(i,j,h))-M*r(i,j,1);
i=8;
E(i,j,h)=(Q2+2/3*F)*(c(i,j+1,h)-c(i,j,h))+beta*Q2*(c(i-1,j,h)-c(i,j,h))+beta*Q2*(c(i+1,j,h)-c(i,j,h))-M*r(i,j,1);
i=9;
E(i,j,h)=Q3*(c(i,j+1,h)-c(i,j,h))+beta*Q2*(c(i-1,j,h)-c(i,j,h))+beta*Q3*(c(i+1,j,h)-c(i,j,h))-M*r(i,j,1);
i=10;
E(i,j,h)=Q3*(c(i,j+1,h)-c(i,j,h))+beta*Q3*(c(i-1,j,h)+c(i,j-1,h)-2*c(i,j,h))-M*r(i,j,1);
i=11;
E(i,j,h)=Q3*c(i-1,j,h)+1/3*F*c(i+1,j,h)-(Q3+1/3*F)*c(i,j,h)+beta*Q3*(c(i,j-1,h)+c(i+1,j,h)-2*c(i,j,h))-M*r(i,j,1);
i=12;
E(i,j,h)=(Q3+F)*(c(i,j-1,h)-c(i,j,h))+beta*Q3*(c(i-1,j,h)-c(i,j,h))+10*beta*Q3*(c(i+1,j,h)-c(i,j,h))-M*r(i,j,1);
i=13;
E(i,j,h)=(Q4-1/3*F)*c(i,j-1,h)+1/3*F*c(i+1,j,h)-Q4*c(i,j,1)+10*beta*Q3*(c(i-1,j,h)-c(i,j,h))+beta*Q4*(c(i+1,j,h)-c(i,j,h))-M*r(i,j,1);
i=14;
E(i,j,h)=(Q4+1/3*F)*(c(i+1,j,h)-c(i,j,h))+beta*Q4*(c(i-1,j,h)+c(i,j-1,h)-2*c(i,j,h))-M*r(i,j,1);
i=15;
E(i,j,h)=(Q4+1/3*F)*(c(i,j+1,h)-c(i,j,h))+beta*Q4*(c(i+1,j,h)+c(i,j-1,h)-2*c(i,j,h))-M*r(i,j,1);
i=16;
E(i,j,h)=(Q4+2/3*F)*(c(i,j+1,h)-c(i,j,h))+beta*Q4*(c(i-1,j,h)-c(i,j,h))-M*r(i,j,1);
%-----------------------------------------------------------------------
j=3;
i=1;
E(i,j,h)=1/3*F*c(i+1,j,h)+(Q1-1/3*F)*c(i,j+1,h)-Q1*c(i,j,h)+beta*Q1*(c(i+1,j,h)-c(i,j,h))-M*r(i,j,1);
i=2;
E(i,j,h)=(Q1+1/3*F)*(c(i+1,j,h)-c(i,j,h))+beta*Q1*(c(i,j+1,h)+c(i-1,j,h)-2*c(i,j,h))-M*r(i,j,1);
i=3;
E(i,j,h)=(Q1+1/3*F)*(c(i,j-1,h)-c(i,j,h))+beta*Q1*(c(i+1,j,h)+c(i,j+1,h)-2*c(i,j,h))-M*r(i,j,1);
i=4;
E(i,j,h)=(Q1+2/3*F)*(c(i,j-1,h)-c(i,j,h))+beta*Q1*(c(i-1,j,h)-c(i,j,h))+10*beta*Q1*(c(i+1,j,h)-c(i,j,h))-M*r(i,j,1);
i=5;
E(i,j,h)=Q2*(c(i,j-1,h)-c(i,j,1))+10*beta*Q1*(c(i-1,j,h)-c(i,j,h))+beta*Q2*(c(i+1,j,h)-c(i,j,h))-M*r(i,j,1);
i=6;
E(i,j,h)=Q2*(c(i,j-1,h)-c(i,j,h))+beta*Q2*(c(i-1,j,h)+c(i,j+1,h)-2*c(i,j,h))-M*r(i,j,1);
i=7;
E(i,j,h)=Q2*c(i-1,j,h)+1/3*F*c(i+1,j,h)-(Q2+1/3*F)*c(i,j,h)+beta*Q2*(c(i,j+1,h)+c(i+1,j,h)-2*c(i,j,h))-M*r(i,j,1);
i=8;
E(i,j,h)=(Q2+F)*(c(i,j+1,h)-c(i,j,h))+beta*Q2*(c(i-1,j,h)-c(i,j,h))+beta*Q2*(c(i+1,j,h)-c(i,j,h))-M*r(i,j,1);
i=9;
E(i,j,h)=(Q3-1/3*F)*c(i,j+1,h)+1/3*F*c(i+1,j,h)-Q3*c(i,j,h)+beta*Q2*(c(i-1,j,h)-c(i,j,h))+beta*Q3*(c(i+1,j,h)-c(i,j,h))-M*r(i,j,1);
i=10;
E(i,j,h)=(Q3+1/3*F)*(c(i+1,j,h)-c(i,j,h))+beta*Q3*(c(i,j+1,h)+c(i-1,j,h)-2*c(i,j,h))-M*r(i,j,1);
i=11;
E(i,j,h)=(Q3+1/3*F)*(c(i,j-1,h)-c(i,j,h))+beta*Q3*(c(i,j+1,h)+c(i+1,j,h)-2*c(i,j,h))-M*r(i,j,1);
i=12;
E(i,j,h)=(Q3+2/3*F)*(c(i,j-1,h)-c(i,j,h))+beta*Q3*(c(i-1,j,h)-c(i,j,h))+10*beta*Q3*(c(i+1,j,h)-c(i,j,h))-M*r(i,j,1);
i=13;
E(i,j,h)=Q4*(c(i,j-1,h)-c(i,j,h))+10*beta*Q3*(c(i-1,j,h)-c(i,j,h))+beta*Q4*(c(i+1,j,h)-c(i,j,h))-M*r(i,j,1);
i=14;
E(i,j,h)=Q4*(c(i,j-1,h)-c(i,j,h))+beta*Q4*(c(i-1,j,h)+c(i,j+1,h)-2*c(i,j,h))-M*r(i,j,1);
i=15;
E(i,j,h)=Q4*c(i-1,j,h)+1/3*F*c(i+1,j,h)-(Q4+1/3*F)*c(i,j,h)+beta*Q4*(c(i,j+1,h)+c(i+1,j,h)-2*c(i,j,h))-M*r(i,j,1);
i=16;
E(i,j,h)=(Q4+F)*(c(i,j+1,h)-c(i,j,h))+beta*Q4*(c(i-1,j,h)-c(i,j,h))-M*r(i,j,1);
%--------------------------------------------------------------------------
j=4;
i=1;
E(i,j,h)=(Q1+2/3*F)*(c(i+1,j,h)-c(i,j,h))-M*r(i,j,1);%出口处平衡方程
i=linspace(2,3,2);
E(i,j,h)=(Q1+2/3*F)*(c(i+1,j,h)-c(i,j,h))+beta*Q1*(c(i,j-1,h)-c(i,j,h))-M*r(i,j,1);
i=4;
E(i,j,h)=(Q1+2/3*F)*(c(i,j-1,h)-c(i,j,h))+10*beta*Q1*(c(i+1,j,h)-c(i,j,h))-M*r(i,j,1);
i=5;
E(i,j,h)=Q2*(c(i,j-1,h)-c(i,j,h))+10*beta*Q1*(c(i-1,j,h)-c(i,j,h))-M*r(i,j,1);
i=linspace(6,7,2);
E(i,j,h)=Q2*(c(i-1,j,h)-c(i,j,h))+beta*Q2*(c(i,j-1,h)-c(i,j,h))-M*r(i,j,1);
i=8;
E(i,j,h)=Q2*c(i-1,j,h)+F*c(i+1,j,h)-(Q2+F)*c(i,j,h)+beta*Q2*(c(i+1,j,h)-c(i,j,h))-M*r(i,j,1);
i=9;
E(i,j,h)=(Q3+2/3*F)*(c(i+1,j,h)-c(i,j,h))+beta*Q2*(c(i-1,j,h)-c(i,j,h))-M*r(i,j,1);
i=linspace(10,11,2);
E(i,j,h)=(Q3+2/3*F)*(c(i+1,j,h)-c(i,j,h))+beta*Q3*(c(i,j-1,h)-c(i,j,h))-M*r(i,j,1);
i=12;
E(i,j,h)=(Q3+2/3*F)*(c(i,j-1,h)-c(i,j,h))+10*beta*Q3*(c(i+1,j,h)-c(i,j,h))-M*r(i,j,1);
i=13;
E(i,j,h)=Q4*(c(i,j-1,h)-c(i,j,h))+10*beta*Q3*(c(i-1,j,h)-c(i,j,h))-M*r(i,j,1);
i=linspace(14,15,2);
E(i,j,h)=Q4*(c(i-1,j,h)-c(i,j,h))+beta*Q4*(c(i,j-1,h)-c(i,j,h))-M*r(i,j,1);
i=16;
E(i,j,h)=F*c0+Q4*c(i-1,j,h)-(F+Q4)*c(i,j,h)-M*r(i,j,1);
%------------------------------------------------------------------
h=2;
j=1;
i=1;
E(i,j,h)=Q1*(c(i,j+1,h)-c(i,j,h))-M*(r(i,j,2)-r(i,j,1));
i=linspace(2,3,2);
E(i,j,h)=Q1*(c(i-1,j,h)-c(i,j,h))+beta*Q1*(c(i,j+1,h)-c(i,j,h))-M*(r(i,j,2)-r(i,j,1));
i=4;
E(i,j,h)=Q1*c(i-1,j,h)+F*c(i+1,j,h)-(Q1+F)*c(i,j,h)+10*beta*Q1*(c(i+1,j,h)-c(i,j,h))-M*(r(i,j,2)-r(i,j,1));
i=5;
E(i,j,h)=(Q2+2/3*F)*(c(i+1,j,h)-c(i,j,h))+10*beta*Q1*(c(i-1,j,h)-c(i,j,h))-M*(r(i,j,2)-r(i,j,1));
i=linspace(6,7,2);
E(i,j,h)=(Q2+2/3*F)*(c(i+1,j,h)-c(i,j,h))+beta*Q2*(c(i,j+1,h)-c(i,j,h))-M*(r(i,j,2)-r(i,j,1));
i=8;
E(i,j,h)=(Q2+2/3*F)*(c(i,j+1,h)-c(i,j,h))+beta*Q2*(c(i+1,j,h)-c(i,j,h))-M*(r(i,j,2)-r(i,j,1));
%第二个循环与第三个循环之间物料交换采用不同的流量以上一个循环为标准
i=9;
E(i,j,h)=Q3*(c(i,j+1,h)-c(i,j,h))+beta*Q2*(c(i-1,j,h)-c(i,j,h))-M*(r(i,j,2)-r(i,j,1));
i=linspace(10,11,2);
E(i,j,h)=Q3*(c(i-1,j,h)-c(i,j,h))+beta*Q3*(c(i,j+1,h)-c(i,j,h))-M*(r(i,j,2)-r(i,j,1));
i=12;
E(i,j,h)=Q3*c(i-1,j,h)+F*c(i+1,j,h)-(Q3+F)*c(i,j,h)+10*beta*Q3*(c(i+1,j,h)-c(i,j,h))-M*(r(i,j,2)-r(i,j,1));
i=13;
E(i,j,h)=(Q4+2/3*F)*(c(i+1,j,h)-c(i,j,h))+10*beta*Q3*(c(i-1,j,h)-c(i,j,h))-M*(r(i,j,2)-r(i,j,1));
i=linspace(14,15,2);
E(i,j,h)=(Q4+2/3*F)*(c(i+1,j,h)-c(i,j,h))+beta*Q4*(c(i,j+1,h)-c(i,j,h))-M*(r(i,j,2)-r(i,j,1));
i=16;
E(i,j,h)=(Q4+2/3*F)*(c(i,j+1,h)-c(i,j,h))-M*(r(i,j,2)-r(i,j,1));
%------------------------------------------------------------------------
j=2;
i=1;
E(i,j,h)=Q1*(c(i,j+1,h)-c(i,j,h))+beta*Q1*(c(i+1,j,h)-c(i,j,h))-M*(r(i,j,2)-r(i,j,1));
i=2;
E(i,j,h)=Q1*(c(i,j+1,h)-c(i,j,h))+beta*Q1*(c(i-1,j,h)+c(i,j-1,h)-2*c(i,j,h))-M*(r(i,j,2)-r(i,j,1));
i=3;
E(i,j,h)=Q1*c(i-1,j,h)+1/3*F*c(i+1,j,h)-(Q1+1/3*F)*c(i,j,h)+beta*Q1*(c(i,j-1,h)+c(i+1,j,h)-2*c(i,j,h))-M*(r(i,j,2)-r(i,j,1));
i=4;
E(i,j,h)=(Q1+F)*(c(i,j-1,h)-c(i,j,h))+10*beta*Q1*(c(i+1,j,h)-c(i,j,h))+beta*Q1*(c(i-1,j,h)-c(i,j,h))-M*(r(i,j,2)-r(i,j,1));
i=5;
E(i,j,h)=(Q2-1/3*F)*c(i,j-1,h)+1/3*F*c(i+1,j,h)-Q2*c(i,j,h)+10*beta*Q1*(c(i-1,j,h)-c(i,j,h))+beta*Q2*(c(i+1,j,h)-c(i,j,h))-M*(r(i,j,2)-r(i,j,1));
i=6;
E(i,j,h)=(Q2+1/3*F)*(c(i+1,j,h)-c(i,j,h))+beta*Q2*(c(i,j-1,h)+c(i-1,j,h)-2*c(i,j,h))-M*(r(i,j,2)-r(i,j,1));
i=7;
E(i,j,h)=(Q2+1/3*F)*(c(i,j+1,h)-c(i,j,h))+beta*Q2*(c(i,j-1,h)+c(i+1,j,h)-2*c(i,j,h))-M*(r(i,j,2)-r(i,j,1));
i=8;
E(i,j,h)=(Q2+2/3*F)*(c(i,j+1,h)-c(i,j,h))+beta*Q2*(c(i-1,j,h)-c(i,j,h))+beta*Q2*(c(i+1,j,h)-c(i,j,h))-M*(r(i,j,2)-r(i,j,1));
i=9;
E(i,j,h)=Q3*(c(i,j+1,h)-c(i,j,h))+beta*Q2*(c(i-1,j,h)-c(i,j,h))+beta*Q3*(c(i+1,j,h)-c(i,j,h))-M*(r(i,j,2)-r(i,j,1));
i=10;
E(i,j,h)=Q3*(c(i,j+1,h)-c(i,j,h))+beta*Q3*(c(i-1,j,h)+c(i,j-1,h)-2*c(i,j,h))-M*(r(i,j,2)-r(i,j,1));
i=11;
E(i,j,h)=Q3*c(i-1,j,h)+1/3*F*c(i+1,j,h)-(Q3+1/3*F)*c(i,j,h)+beta*Q3*(c(i,j-1,h)+c(i+1,j,h)-2*c(i,j,h))-M*(r(i,j,2)-r(i,j,1));
i=12;
E(i,j,h)=(Q3+F)*(c(i,j-1,h)-c(i,j,h))+beta*Q3*(c(i-1,j,h)-c(i,j,h))+10*beta*Q3*(c(i+1,j,h)-c(i,j,h))-M*(r(i,j,2)-r(i,j,1));
i=13;
E(i,j,h)=(Q4-1/3*F)*c(i,j-1,h)+1/3*F*c(i+1,j,h)-Q4*c(i,j,1)+10*beta*Q3*(c(i-1,j,h)-c(i,j,h))+beta*Q4*(c(i+1,j,h)-c(i,j,h))-M*(r(i,j,2)-r(i,j,1));
i=14;
E(i,j,h)=(Q4+1/3*F)*(c(i+1,j,h)-c(i,j,h))+beta*Q4*(c(i-1,j,h)+c(i,j-1,h)-2*c(i,j,h))-M*(r(i,j,2)-r(i,j,1));
i=15;
E(i,j,h)=(Q4+1/3*F)*(c(i,j+1,h)-c(i,j,h))+beta*Q4*(c(i+1,j,h)+c(i,j-1,h)-2*c(i,j,h))-M*(r(i,j,2)-r(i,j,1));
i=16;
E(i,j,h)=(Q4+2/3*F)*(c(i,j+1,h)-c(i,j,h))+beta*Q4*(c(i-1,j,h)-c(i,j,h))-M*(r(i,j,2)-r(i,j,1));
%-----------------------------------------------------------------------
j=3;
i=1;
E(i,j,h)=1/3*F*c(i+1,j,h)+(Q1-1/3*F)*c(i,j+1,h)-Q1*c(i,j,h)+beta*Q1*(c(i+1,j,h)-c(i,j,h))-M*(r(i,j,2)-r(i,j,1));
i=2;
E(i,j,h)=(Q1+1/3*F)*(c(i+1,j,h)-c(i,j,h))+beta*Q1*(c(i,j+1,h)+c(i-1,j,h)-2*c(i,j,h))-M*(r(i,j,2)-r(i,j,1));
i=3;
E(i,j,h)=(Q1+1/3*F)*(c(i,j-1,h)-c(i,j,h))+beta*Q1*(c(i+1,j,h)+c(i,j+1,h)-2*c(i,j,h))-M*(r(i,j,2)-r(i,j,1));
i=4;
E(i,j,h)=(Q1+2/3*F)*(c(i,j-1,h)-c(i,j,h))+beta*Q1*(c(i-1,j,h)-c(i,j,h))+10*beta*Q1*(c(i+1,j,h)-c(i,j,h))-M*(r(i,j,2)-r(i,j,1));
i=5;
E(i,j,h)=Q2*(c(i,j-1,h)-c(i,j,1))+10*beta*Q1*(c(i-1,j,h)-c(i,j,h))+beta*Q2*(c(i+1,j,h)-c(i,j,h))-M*(r(i,j,2)-r(i,j,1));
i=6;
E(i,j,h)=Q2*(c(i,j-1,h)-c(i,j,h))+beta*Q2*(c(i-1,j,h)+c(i,j+1,h)-2*c(i,j,h))-M*(r(i,j,2)-r(i,j,1));
i=7;
E(i,j,h)=Q2*c(i-1,j,h)+1/3*F*c(i+1,j,h)-(Q2+1/3*F)*c(i,j,h)+beta*Q2*(c(i,j+1,h)+c(i+1,j,h)-2*c(i,j,h))-M*(r(i,j,2)-r(i,j,1));
i=8;
E(i,j,h)=(Q2+F)*(c(i,j+1,h)-c(i,j,h))+beta*Q2*(c(i-1,j,h)-c(i,j,h))+beta*Q2*(c(i+1,j,h)-c(i,j,h))-M*(r(i,j,2)-r(i,j,1));
i=9;
E(i,j,h)=(Q3-1/3*F)*c(i,j+1,h)+1/3*F*c(i+1,j,h)-Q3*c(i,j,h)+beta*Q2*(c(i-1,j,h)-c(i,j,h))+beta*Q3*(c(i+1,j,h)-c(i,j,h))-M*(r(i,j,2)-r(i,j,1));
i=10;
E(i,j,h)=(Q3+1/3*F)*(c(i+1,j,h)-c(i,j,h))+beta*Q3*(c(i,j+1,h)+c(i-1,j,h)-2*c(i,j,h))-M*(r(i,j,2)-r(i,j,1));
i=11;
E(i,j,h)=(Q3+1/3*F)*(c(i,j-1,h)-c(i,j,h))+beta*Q3*(c(i,j+1,h)+c(i+1,j,h)-2*c(i,j,h))-M*(r(i,j,2)-r(i,j,1));
i=12;
E(i,j,h)=(Q3+2/3*F)*(c(i,j-1,h)-c(i,j,h))+beta*Q3*(c(i-1,j,h)-c(i,j,h))+10*beta*Q3*(c(i+1,j,h)-c(i,j,h))-M*(r(i,j,2)-r(i,j,1));
i=13;
E(i,j,h)=Q4*(c(i,j-1,h)-c(i,j,h))+10*beta*Q3*(c(i-1,j,h)-c(i,j,h))+beta*Q4*(c(i+1,j,h)-c(i,j,h))-M*(r(i,j,2)-r(i,j,1));
i=14;
E(i,j,h)=Q4*(c(i,j-1,h)-c(i,j,h))+beta*Q4*(c(i-1,j,h)+c(i,j+1,h)-2*c(i,j,h))-M*(r(i,j,2)-r(i,j,1));
i=15;
E(i,j,h)=Q4*c(i-1,j,h)+1/3*F*c(i+1,j,h)-(Q4+1/3*F)*c(i,j,h)+beta*Q4*(c(i,j+1,h)+c(i+1,j,h)-2*c(i,j,h))-M*(r(i,j,2)-r(i,j,1));
i=16;
E(i,j,h)=(Q4+F)*(c(i,j+1,h)-c(i,j,h))+beta*Q4*(c(i-1,j,h)-c(i,j,h))-M*(r(i,j,2)-r(i,j,1));
%--------------------------------------------------------------------------
j=4;
i=1;
E(i,j,h)=(Q1+2/3*F)*(c(i+1,j,h)-c(i,j,h))-M*(r(i,j,2)-r(i,j,1));%出口处平衡方程
i=linspace(2,3,2);
E(i,j,h)=(Q1+2/3*F)*(c(i+1,j,h)-c(i,j,h))+beta*Q1*(c(i,j-1,h)-c(i,j,h))-M*(r(i,j,2)-r(i,j,1));
i=4;
E(i,j,h)=(Q1+2/3*F)*(c(i,j-1,h)-c(i,j,h))+10*beta*Q1*(c(i+1,j,h)-c(i,j,h))-M*(r(i,j,2)-r(i,j,1));
i=5;
E(i,j,h)=Q2*(c(i,j-1,h)-c(i,j,h))+10*beta*Q1*(c(i-1,j,h)-c(i,j,h))-M*(r(i,j,2)-r(i,j,1));
i=linspace(6,7,2);
E(i,j,h)=Q2*(c(i-1,j,h)-c(i,j,h))+beta*Q2*(c(i,j-1,h)-c(i,j,h))-M*(r(i,j,2)-r(i,j,1));
i=8;
E(i,j,h)=Q2*c(i-1,j,h)+F*c(i+1,j,h)-(Q2+F)*c(i,j,h)+beta*Q2*(c(i+1,j,h)-c(i,j,h))-M*(r(i,j,2)-r(i,j,1));
i=9;
E(i,j,h)=(Q3+2/3*F)*(c(i+1,j,h)-c(i,j,h))+beta*Q2*(c(i-1,j,h)-c(i,j,h))-M*(r(i,j,2)-r(i,j,1));
i=linspace(10,11,2);
E(i,j,h)=(Q3+2/3*F)*(c(i+1,j,h)-c(i,j,h))+beta*Q3*(c(i,j-1,h)-c(i,j,h))-M*(r(i,j,2)-r(i,j,1));
i=12;
E(i,j,h)=(Q3+2/3*F)*(c(i,j-1,h)-c(i,j,h))+10*beta*Q3*(c(i+1,j,h)-c(i,j,h))-M*(r(i,j,2)-r(i,j,1));
i=13;
E(i,j,h)=Q4*(c(i,j-1,h)-c(i,j,h))+10*beta*Q3*(c(i-1,j,h)-c(i,j,h))-M*(r(i,j,2)-r(i,j,1));
i=linspace(14,15,2);
E(i,j,h)=Q4*(c(i-1,j,h)-c(i,j,h))+beta*Q4*(c(i,j-1,h)-c(i,j,h))-M*(r(i,j,2)-r(i,j,1));
i=16;
E(i,j,h)=F*0+Q4*c(i-1,j,h)-(F+Q4)*c(i,j,h)-M*(r(i,j,2)-r(i,j,1));
%------------------------------------------------------------------
h=3;
j=1;
i=1;
E(i,j,h)=Q1*(c(i,j+1,h)-c(i,j,h))-M*(r(i,j,3)-r(i,j,2));
i=linspace(2,3,2);
E(i,j,h)=Q1*(c(i-1,j,h)-c(i,j,h))+beta*Q1*(c(i,j+1,h)-c(i,j,h))-M*(r(i,j,3)-r(i,j,2));
i=4;
E(i,j,h)=Q1*c(i-1,j,h)+F*c(i+1,j,h)-(Q1+F)*c(i,j,h)+10*beta*Q1*(c(i+1,j,h)-c(i,j,h))-M*(r(i,j,3)-r(i,j,2));
i=5;
E(i,j,h)=(Q2+2/3*F)*(c(i+1,j,h)-c(i,j,h))+10*beta*Q1*(c(i-1,j,h)-c(i,j,h))-M*(r(i,j,3)-r(i,j,2));
i=linspace(6,7,2);
E(i,j,h)=(Q2+2/3*F)*(c(i+1,j,h)-c(i,j,h))+beta*Q2*(c(i,j+1,h)-c(i,j,h))-M*(r(i,j,3)-r(i,j,2));
i=8;
E(i,j,h)=(Q2+2/3*F)*(c(i,j+1,h)-c(i,j,h))+beta*Q2*(c(i+1,j,h)-c(i,j,h))-M*(r(i,j,3)-r(i,j,2));
%第二个循环与第三个循环之间物料交换采用不同的流量以上一个循环为标准
i=9;
E(i,j,h)=Q3*(c(i,j+1,h)-c(i,j,h))+beta*Q2*(c(i-1,j,h)-c(i,j,h))-M*(r(i,j,3)-r(i,j,2));
i=linspace(10,11,2);
E(i,j,h)=Q3*(c(i-1,j,h)-c(i,j,h))+beta*Q3*(c(i,j+1,h)-c(i,j,h))-M*(r(i,j,3)-r(i,j,2));
i=12;
E(i,j,h)=Q3*c(i-1,j,h)+F*c(i+1,j,h)-(Q3+F)*c(i,j,h)+10*beta*Q3*(c(i+1,j,h)-c(i,j,h))-M*(r(i,j,3)-r(i,j,2));
i=13;
E(i,j,h)=(Q4+2/3*F)*(c(i+1,j,h)-c(i,j,h))+10*beta*Q3*(c(i-1,j,h)-c(i,j,h))-M*(r(i,j,3)-r(i,j,2));
i=linspace(14,15,2);
E(i,j,h)=(Q4+2/3*F)*(c(i+1,j,h)-c(i,j,h))+beta*Q4*(c(i,j+1,h)-c(i,j,h))-M*(r(i,j,3)-r(i,j,2));
i=16;
E(i,j,h)=(Q4+2/3*F)*(c(i,j+1,h)-c(i,j,h))-M*(r(i,j,3)-r(i,j,2));
%------------------------------------------------------------------------
j=2;
i=1;
E(i,j,h)=Q1*(c(i,j+1,h)-c(i,j,h))+beta*Q1*(c(i+1,j,h)-c(i,j,h))-M*(r(i,j,3)-r(i,j,2));
i=2;
E(i,j,h)=Q1*(c(i,j+1,h)-c(i,j,h))+beta*Q1*(c(i-1,j,h)+c(i,j-1,h)-2*c(i,j,h))-M*(r(i,j,3)-r(i,j,2));
i=3;
E(i,j,h)=Q1*c(i-1,j,h)+1/3*F*c(i+1,j,h)-(Q1+1/3*F)*c(i,j,h)+beta*Q1*(c(i,j-1,h)+c(i+1,j,h)-2*c(i,j,h))-M*(r(i,j,3)-r(i,j,2));
i=4;
E(i,j,h)=(Q1+F)*(c(i,j-1,h)-c(i,j,h))+10*beta*Q1*(c(i+1,j,h)-c(i,j,h))+beta*Q1*(c(i-1,j,h)-c(i,j,h))-M*(r(i,j,3)-r(i,j,2));
i=5;
E(i,j,h)=(Q2-1/3*F)*c(i,j-1,h)+1/3*F*c(i+1,j,h)-Q2*c(i,j,h)+10*beta*Q1*(c(i-1,j,h)-c(i,j,h))+beta*Q2*(c(i+1,j,h)-c(i,j,h))-M*(r(i,j,3)-r(i,j,2));
i=6;
E(i,j,h)=(Q2+1/3*F)*(c(i+1,j,h)-c(i,j,h))+beta*Q2*(c(i,j-1,h)+c(i-1,j,h)-2*c(i,j,h))-M*(r(i,j,3)-r(i,j,2));
i=7;
E(i,j,h)=(Q2+1/3*F)*(c(i,j+1,h)-c(i,j,h))+beta*Q2*(c(i,j-1,h)+c(i+1,j,h)-2*c(i,j,h))-M*(r(i,j,3)-r(i,j,2));
i=8;
E(i,j,h)=(Q2+2/3*F)*(c(i,j+1,h)-c(i,j,h))+beta*Q2*(c(i-1,j,h)-c(i,j,h))+beta*Q2*(c(i+1,j,h)-c(i,j,h))-M*(r(i,j,3)-r(i,j,2));
i=9;
E(i,j,h)=Q3*(c(i,j+1,h)-c(i,j,h))+beta*Q2*(c(i-1,j,h)-c(i,j,h))+beta*Q3*(c(i+1,j,h)-c(i,j,h))-M*(r(i,j,3)-r(i,j,2));
i=10;
E(i,j,h)=Q3*(c(i,j+1,h)-c(i,j,h))+beta*Q3*(c(i-1,j,h)+c(i,j-1,h)-2*c(i,j,h))-M*(r(i,j,3)-r(i,j,2));
i=11;
E(i,j,h)=Q3*c(i-1,j,h)+1/3*F*c(i+1,j,h)-(Q3+1/3*F)*c(i,j,h)+beta*Q3*(c(i,j-1,h)+c(i+1,j,h)-2*c(i,j,h))-M*(r(i,j,3)-r(i,j,2));
i=12;
E(i,j,h)=(Q3+F)*(c(i,j-1,h)-c(i,j,h))+beta*Q3*(c(i-1,j,h)-c(i,j,h))+10*beta*Q3*(c(i+1,j,h)-c(i,j,h))-M*(r(i,j,3)-r(i,j,2));
i=13;
E(i,j,h)=(Q4-1/3*F)*c(i,j-1,h)+1/3*F*c(i+1,j,h)-Q4*c(i,j,1)+10*beta*Q3*(c(i-1,j,h)-c(i,j,h))+beta*Q4*(c(i+1,j,h)-c(i,j,h))-M*(r(i,j,3)-r(i,j,2));
i=14;
E(i,j,h)=(Q4+1/3*F)*(c(i+1,j,h)-c(i,j,h))+beta*Q4*(c(i-1,j,h)+c(i,j-1,h)-2*c(i,j,h))-M*(r(i,j,3)-r(i,j,2));
i=15;
E(i,j,h)=(Q4+1/3*F)*(c(i,j+1,h)-c(i,j,h))+beta*Q4*(c(i+1,j,h)+c(i,j-1,h)-2*c(i,j,h))-M*(r(i,j,3)-r(i,j,2));
i=16;
E(i,j,h)=(Q4+2/3*F)*(c(i,j+1,h)-c(i,j,h))+beta*Q4*(c(i-1,j,h)-c(i,j,h))-M*(r(i,j,3)-r(i,j,2));
%-----------------------------------------------------------------------
j=3;
i=1;
E(i,j,h)=1/3*F*c(i+1,j,h)+(Q1-1/3*F)*c(i,j+1,h)-Q1*c(i,j,h)+beta*Q1*(c(i+1,j,h)-c(i,j,h))-M*(r(i,j,3)-r(i,j,2));
i=2;
E(i,j,h)=(Q1+1/3*F)*(c(i+1,j,h)-c(i,j,h))+beta*Q1*(c(i,j+1,h)+c(i-1,j,h)-2*c(i,j,h))-M*(r(i,j,3)-r(i,j,2));
i=3;
E(i,j,h)=(Q1+1/3*F)*(c(i,j-1,h)-c(i,j,h))+beta*Q1*(c(i+1,j,h)+c(i,j+1,h)-2*c(i,j,h))-M*(r(i,j,3)-r(i,j,2));
i=4;
E(i,j,h)=(Q1+2/3*F)*(c(i,j-1,h)-c(i,j,h))+beta*Q1*(c(i-1,j,h)-c(i,j,h))+10*beta*Q1*(c(i+1,j,h)-c(i,j,h))-M*(r(i,j,3)-r(i,j,2));
i=5;
E(i,j,h)=Q2*(c(i,j-1,h)-c(i,j,1))+10*beta*Q1*(c(i-1,j,h)-c(i,j,h))+beta*Q2*(c(i+1,j,h)-c(i,j,h))-M*(r(i,j,3)-r(i,j,2));
i=6;
E(i,j,h)=Q2*(c(i,j-1,h)-c(i,j,h))+beta*Q2*(c(i-1,j,h)+c(i,j+1,h)-2*c(i,j,h))-M*(r(i,j,3)-r(i,j,2));
i=7;
E(i,j,h)=Q2*c(i-1,j,h)+1/3*F*c(i+1,j,h)-(Q2+1/3*F)*c(i,j,h)+beta*Q2*(c(i,j+1,h)+c(i+1,j,h)-2*c(i,j,h))-M*(r(i,j,3)-r(i,j,2));
i=8;
E(i,j,h)=(Q2+F)*(c(i,j+1,h)-c(i,j,h))+beta*Q2*(c(i-1,j,h)-c(i,j,h))+beta*Q2*(c(i+1,j,h)-c(i,j,h))-M*(r(i,j,3)-r(i,j,2));
i=9;
E(i,j,h)=(Q3-1/3*F)*c(i,j+1,h)+1/3*F*c(i+1,j,h)-Q3*c(i,j,h)+beta*Q2*(c(i-1,j,h)-c(i,j,h))+beta*Q3*(c(i+1,j,h)-c(i,j,h))-M*(r(i,j,3)-r(i,j,2));
i=10;
E(i,j,h)=(Q3+1/3*F)*(c(i+1,j,h)-c(i,j,h))+beta*Q3*(c(i,j+1,h)+c(i-1,j,h)-2*c(i,j,h))-M*(r(i,j,3)-r(i,j,2));
i=11;
E(i,j,h)=(Q3+1/3*F)*(c(i,j-1,h)-c(i,j,h))+beta*Q3*(c(i,j+1,h)+c(i+1,j,h)-2*c(i,j,h))-M*(r(i,j,3)-r(i,j,2));
i=12;
E(i,j,h)=(Q3+2/3*F)*(c(i,j-1,h)-c(i,j,h))+beta*Q3*(c(i-1,j,h)-c(i,j,h))+10*beta*Q3*(c(i+1,j,h)-c(i,j,h))-M*(r(i,j,3)-r(i,j,2));
i=13;
E(i,j,h)=Q4*(c(i,j-1,h)-c(i,j,h))+10*beta*Q3*(c(i-1,j,h)-c(i,j,h))+beta*Q4*(c(i+1,j,h)-c(i,j,h))-M*(r(i,j,3)-r(i,j,2));
i=14;
E(i,j,h)=Q4*(c(i,j-1,h)-c(i,j,h))+beta*Q4*(c(i-1,j,h)+c(i,j+1,h)-2*c(i,j,h))-M*(r(i,j,3)-r(i,j,2));
i=15;
E(i,j,h)=Q4*c(i-1,j,h)+1/3*F*c(i+1,j,h)-(Q4+1/3*F)*c(i,j,h)+beta*Q4*(c(i,j+1,h)+c(i+1,j,h)-2*c(i,j,h))-M*(r(i,j,3)-r(i,j,2));
i=16;
E(i,j,h)=(Q4+F)*(c(i,j+1,h)-c(i,j,h))+beta*Q4*(c(i-1,j,h)-c(i,j,h))-M*(r(i,j,3)-r(i,j,2));
%--------------------------------------------------------------------------
j=4;
i=1;
E(i,j,h)=(Q1+2/3*F)*(c(i+1,j,h)-c(i,j,h))-M*(r(i,j,3)-r(i,j,2));%出口处平衡方程
i=linspace(2,3,2);
E(i,j,h)=(Q1+2/3*F)*(c(i+1,j,h)-c(i,j,h))+beta*Q1*(c(i,j-1,h)-c(i,j,h))-M*(r(i,j,3)-r(i,j,2));
i=4;
E(i,j,h)=(Q1+2/3*F)*(c(i,j-1,h)-c(i,j,h))+10*beta*Q1*(c(i+1,j,h)-c(i,j,h))-M*(r(i,j,3)-r(i,j,2));
i=5;
E(i,j,h)=Q2*(c(i,j-1,h)-c(i,j,h))+10*beta*Q1*(c(i-1,j,h)-c(i,j,h))-M*(r(i,j,3)-r(i,j,2));
i=linspace(6,7,2);
E(i,j,h)=Q2*(c(i-1,j,h)-c(i,j,h))+beta*Q2*(c(i,j-1,h)-c(i,j,h))-M*(r(i,j,3)-r(i,j,2));
i=8;
E(i,j,h)=Q2*c(i-1,j,h)+F*c(i+1,j,h)-(Q2+F)*c(i,j,h)+beta*Q2*(c(i+1,j,h)-c(i,j,h))-M*(r(i,j,3)-r(i,j,2));
i=9;
E(i,j,h)=(Q3+2/3*F)*(c(i+1,j,h)-c(i,j,h))+beta*Q2*(c(i-1,j,h)-c(i,j,h))-M*(r(i,j,3)-r(i,j,2));
i=linspace(10,11,2);
E(i,j,h)=(Q3+2/3*F)*(c(i+1,j,h)-c(i,j,h))+beta*Q3*(c(i,j-1,h)-c(i,j,h))-M*(r(i,j,3)-r(i,j,2));
i=12;
E(i,j,h)=(Q3+2/3*F)*(c(i,j-1,h)-c(i,j,h))+10*beta*Q3*(c(i+1,j,h)-c(i,j,h))-M*(r(i,j,3)-r(i,j,2));
i=13;
E(i,j,h)=Q4*(c(i,j-1,h)-c(i,j,h))+10*beta*Q3*(c(i-1,j,h)-c(i,j,h))-M*(r(i,j,3)-r(i,j,2));
i=linspace(14,15,2);
E(i,j,h)=Q4*(c(i-1,j,h)-c(i,j,h))+beta*Q4*(c(i,j-1,h)-c(i,j,h))-M*(r(i,j,3)-r(i,j,2));
i=16;
E(i,j,h)=F*0+Q4*c(i-1,j,h)-(F+Q4)*c(i,j,h)-M*(r(i,j,3)-r(i,j,2));
%------------------------------------------------------------------
h=4;
j=1;
i=1;
E(i,j,h)=Q1*(c(i,j+1,h)-c(i,j,h))-M*(r(i,j,4)-r(i,j,3));
i=linspace(2,3,2);
E(i,j,h)=Q1*(c(i-1,j,h)-c(i,j,h))+beta*Q1*(c(i,j+1,h)-c(i,j,h))-M*(r(i,j,4)-r(i,j,3));
i=4;
E(i,j,h)=Q1*c(i-1,j,h)+F*c(i+1,j,h)-(Q1+F)*c(i,j,h)+10*beta*Q1*(c(i+1,j,h)-c(i,j,h))-M*(r(i,j,4)-r(i,j,3));
i=5;
E(i,j,h)=(Q2+2/3*F)*(c(i+1,j,h)-c(i,j,h))+10*beta*Q1*(c(i-1,j,h)-c(i,j,h))-M*(r(i,j,4)-r(i,j,3));
i=linspace(6,7,2);
E(i,j,h)=(Q2+2/3*F)*(c(i+1,j,h)-c(i,j,h))+beta*Q2*(c(i,j+1,h)-c(i,j,h))-M*(r(i,j,4)-r(i,j,3));
i=8;
E(i,j,h)=(Q2+2/3*F)*(c(i,j+1,h)-c(i,j,h))+beta*Q2*(c(i+1,j,h)-c(i,j,h))-M*(r(i,j,4)-r(i,j,3));
%第二个循环与第三个循环之间物料交换采用不同的流量以上一个循环为标准
i=9;
E(i,j,h)=Q3*(c(i,j+1,h)-c(i,j,h))+beta*Q2*(c(i-1,j,h)-c(i,j,h))-M*(r(i,j,4)-r(i,j,3));
i=linspace(10,11,2);
E(i,j,h)=Q3*(c(i-1,j,h)-c(i,j,h))+beta*Q3*(c(i,j+1,h)-c(i,j,h))-M*(r(i,j,4)-r(i,j,3));
i=12;
E(i,j,h)=Q3*c(i-1,j,h)+F*c(i+1,j,h)-(Q3+F)*c(i,j,h)+10*beta*Q3*(c(i+1,j,h)-c(i,j,h))-M*(r(i,j,4)-r(i,j,3));
i=13;
E(i,j,h)=(Q4+2/3*F)*(c(i+1,j,h)-c(i,j,h))+10*beta*Q3*(c(i-1,j,h)-c(i,j,h))-M*(r(i,j,4)-r(i,j,3));
i=linspace(14,15,2);
E(i,j,h)=(Q4+2/3*F)*(c(i+1,j,h)-c(i,j,h))+beta*Q4*(c(i,j+1,h)-c(i,j,h))-M*(r(i,j,4)-r(i,j,3));
i=16;
E(i,j,h)=(Q4+2/3*F)*(c(i,j+1,h)-c(i,j,h))-M*(r(i,j,4)-r(i,j,3));
%------------------------------------------------------------------------
j=2;
i=1;
E(i,j,h)=Q1*(c(i,j+1,h)-c(i,j,h))+beta*Q1*(c(i+1,j,h)-c(i,j,h))-M*(r(i,j,4)-r(i,j,3));
i=2;
E(i,j,h)=Q1*(c(i,j+1,h)-c(i,j,h))+beta*Q1*(c(i-1,j,h)+c(i,j-1,h)-2*c(i,j,h))-M*(r(i,j,4)-r(i,j,3));
i=3;
E(i,j,h)=Q1*c(i-1,j,h)+1/3*F*c(i+1,j,h)-(Q1+1/3*F)*c(i,j,h)+beta*Q1*(c(i,j-1,h)+c(i+1,j,h)-2*c(i,j,h))-M*(r(i,j,4)-r(i,j,3));
i=4;
E(i,j,h)=(Q1+F)*(c(i,j-1,h)-c(i,j,h))+10*beta*Q1*(c(i+1,j,h)-c(i,j,h))+beta*Q1*(c(i-1,j,h)-c(i,j,h))-M*(r(i,j,4)-r(i,j,3));
i=5;
E(i,j,h)=(Q2-1/3*F)*c(i,j-1,h)+1/3*F*c(i+1,j,h)-Q2*c(i,j,h)+10*beta*Q1*(c(i-1,j,h)-c(i,j,h))+beta*Q2*(c(i+1,j,h)-c(i,j,h))-M*(r(i,j,4)-r(i,j,3));
i=6;
E(i,j,h)=(Q2+1/3*F)*(c(i+1,j,h)-c(i,j,h))+beta*Q2*(c(i,j-1,h)+c(i-1,j,h)-2*c(i,j,h))-M*(r(i,j,4)-r(i,j,3));
i=7;
E(i,j,h)=(Q2+1/3*F)*(c(i,j+1,h)-c(i,j,h))+beta*Q2*(c(i,j-1,h)+c(i+1,j,h)-2*c(i,j,h))-M*(r(i,j,4)-r(i,j,3));
i=8;
E(i,j,h)=(Q2+2/3*F)*(c(i,j+1,h)-c(i,j,h))+beta*Q2*(c(i-1,j,h)-c(i,j,h))+beta*Q2*(c(i+1,j,h)-c(i,j,h))-M*(r(i,j,4)-r(i,j,3));
i=9;
E(i,j,h)=Q3*(c(i,j+1,h)-c(i,j,h))+beta*Q2*(c(i-1,j,h)-c(i,j,h))+beta*Q3*(c(i+1,j,h)-c(i,j,h))-M*(r(i,j,4)-r(i,j,3));
i=10;
E(i,j,h)=Q3*(c(i,j+1,h)-c(i,j,h))+beta*Q3*(c(i-1,j,h)+c(i,j-1,h)-2*c(i,j,h))-M*(r(i,j,4)-r(i,j,3));
i=11;
E(i,j,h)=Q3*c(i-1,j,h)+1/3*F*c(i+1,j,h)-(Q3+1/3*F)*c(i,j,h)+beta*Q3*(c(i,j-1,h)+c(i+1,j,h)-2*c(i,j,h))-M*(r(i,j,4)-r(i,j,3));
i=12;
E(i,j,h)=(Q3+F)*(c(i,j-1,h)-c(i,j,h))+beta*Q3*(c(i-1,j,h)-c(i,j,h))+10*beta*Q3*(c(i+1,j,h)-c(i,j,h))-M*(r(i,j,4)-r(i,j,3));
i=13;
E(i,j,h)=(Q4-1/3*F)*c(i,j-1,h)+1/3*F*c(i+1,j,h)-Q4*c(i,j,1)+10*beta*Q3*(c(i-1,j,h)-c(i,j,h))+beta*Q4*(c(i+1,j,h)-c(i,j,h))-M*(r(i,j,4)-r(i,j,3));
i=14;
E(i,j,h)=(Q4+1/3*F)*(c(i+1,j,h)-c(i,j,h))+beta*Q4*(c(i-1,j,h)+c(i,j-1,h)-2*c(i,j,h))-M*(r(i,j,4)-r(i,j,3));
i=15;
E(i,j,h)=(Q4+1/3*F)*(c(i,j+1,h)-c(i,j,h))+beta*Q4*(c(i+1,j,h)+c(i,j-1,h)-2*c(i,j,h))-M*(r(i,j,4)-r(i,j,3));
i=16;
E(i,j,h)=(Q4+2/3*F)*(c(i,j+1,h)-c(i,j,h))+beta*Q4*(c(i-1,j,h)-c(i,j,h))-M*(r(i,j,4)-r(i,j,3));
%-----------------------------------------------------------------------
j=3;
i=1;
E(i,j,h)=1/3*F*c(i+1,j,h)+(Q1-1/3*F)*c(i,j+1,h)-Q1*c(i,j,h)+beta*Q1*(c(i+1,j,h)-c(i,j,h))-M*(r(i,j,4)-r(i,j,3));
i=2;
E(i,j,h)=(Q1+1/3*F)*(c(i+1,j,h)-c(i,j,h))+beta*Q1*(c(i,j+1,h)+c(i-1,j,h)-2*c(i,j,h))-M*(r(i,j,4)-r(i,j,3));
i=3;
E(i,j,h)=(Q1+1/3*F)*(c(i,j-1,h)-c(i,j,h))+beta*Q1*(c(i+1,j,h)+c(i,j+1,h)-2*c(i,j,h))-M*(r(i,j,4)-r(i,j,3));
i=4;
E(i,j,h)=(Q1+2/3*F)*(c(i,j-1,h)-c(i,j,h))+beta*Q1*(c(i-1,j,h)-c(i,j,h))+10*beta*Q1*(c(i+1,j,h)-c(i,j,h))-M*(r(i,j,4)-r(i,j,3));
i=5;
E(i,j,h)=Q2*(c(i,j-1,h)-c(i,j,1))+10*beta*Q1*(c(i-1,j,h)-c(i,j,h))+beta*Q2*(c(i+1,j,h)-c(i,j,h))-M*(r(i,j,4)-r(i,j,3));
i=6;
E(i,j,h)=Q2*(c(i,j-1,h)-c(i,j,h))+beta*Q2*(c(i-1,j,h)+c(i,j+1,h)-2*c(i,j,h))-M*(r(i,j,4)-r(i,j,3));
i=7;
E(i,j,h)=Q2*c(i-1,j,h)+1/3*F*c(i+1,j,h)-(Q2+1/3*F)*c(i,j,h)+beta*Q2*(c(i,j+1,h)+c(i+1,j,h)-2*c(i,j,h))-M*(r(i,j,4)-r(i,j,3));
i=8;
E(i,j,h)=(Q2+F)*(c(i,j+1,h)-c(i,j,h))+beta*Q2*(c(i-1,j,h)-c(i,j,h))+beta*Q2*(c(i+1,j,h)-c(i,j,h))-M*(r(i,j,4)-r(i,j,3));
i=9;
E(i,j,h)=(Q3-1/3*F)*c(i,j+1,h)+1/3*F*c(i+1,j,h)-Q3*c(i,j,h)+beta*Q2*(c(i-1,j,h)-c(i,j,h))+beta*Q3*(c(i+1,j,h)-c(i,j,h))-M*(r(i,j,4)-r(i,j,3));
i=10;
E(i,j,h)=(Q3+1/3*F)*(c(i+1,j,h)-c(i,j,h))+beta*Q3*(c(i,j+1,h)+c(i-1,j,h)-2*c(i,j,h))-M*(r(i,j,4)-r(i,j,3));
i=11;
E(i,j,h)=(Q3+1/3*F)*(c(i,j-1,h)-c(i,j,h))+beta*Q3*(c(i,j+1,h)+c(i+1,j,h)-2*c(i,j,h))-M*(r(i,j,4)-r(i,j,3));
i=12;
E(i,j,h)=(Q3+2/3*F)*(c(i,j-1,h)-c(i,j,h))+beta*Q3*(c(i-1,j,h)-c(i,j,h))+10*beta*Q3*(c(i+1,j,h)-c(i,j,h))-M*(r(i,j,4)-r(i,j,3));
i=13;
E(i,j,h)=Q4*(c(i,j-1,h)-c(i,j,h))+10*beta*Q3*(c(i-1,j,h)-c(i,j,h))+beta*Q4*(c(i+1,j,h)-c(i,j,h))-M*(r(i,j,4)-r(i,j,3));
i=14;
E(i,j,h)=Q4*(c(i,j-1,h)-c(i,j,h))+beta*Q4*(c(i-1,j,h)+c(i,j+1,h)-2*c(i,j,h))-M*(r(i,j,4)-r(i,j,3));
i=15;
E(i,j,h)=Q4*c(i-1,j,h)+1/3*F*c(i+1,j,h)-(Q4+1/3*F)*c(i,j,h)+beta*Q4*(c(i,j+1,h)+c(i+1,j,h)-2*c(i,j,h))-M*(r(i,j,4)-r(i,j,3));
i=16;
E(i,j,h)=(Q4+F)*(c(i,j+1,h)-c(i,j,h))+beta*Q4*(c(i-1,j,h)-c(i,j,h))-M*(r(i,j,4)-r(i,j,3));
%--------------------------------------------------------------------------
j=4;
i=1;
E(i,j,h)=(Q1+2/3*F)*(c(i+1,j,h)-c(i,j,h))-M*(r(i,j,4)-r(i,j,3));%出口处平衡方程
i=linspace(2,3,2);
E(i,j,h)=(Q1+2/3*F)*(c(i+1,j,h)-c(i,j,h))+beta*Q1*(c(i,j-1,h)-c(i,j,h))-M*(r(i,j,4)-r(i,j,3));
i=4;
E(i,j,h)=(Q1+2/3*F)*(c(i,j-1,h)-c(i,j,h))+10*beta*Q1*(c(i+1,j,h)-c(i,j,h))-M*(r(i,j,4)-r(i,j,3));
i=5;
E(i,j,h)=Q2*(c(i,j-1,h)-c(i,j,h))+10*beta*Q1*(c(i-1,j,h)-c(i,j,h))-M*(r(i,j,4)-r(i,j,3));
i=linspace(6,7,2);
E(i,j,h)=Q2*(c(i-1,j,h)-c(i,j,h))+beta*Q2*(c(i,j-1,h)-c(i,j,h))-M*(r(i,j,4)-r(i,j,3));
i=8;
E(i,j,h)=Q2*c(i-1,j,h)+F*c(i+1,j,h)-(Q2+F)*c(i,j,h)+beta*Q2*(c(i+1,j,h)-c(i,j,h))-M*(r(i,j,4)-r(i,j,3));
i=9;
E(i,j,h)=(Q3+2/3*F)*(c(i+1,j,h)-c(i,j,h))+beta*Q2*(c(i-1,j,h)-c(i,j,h))-M*(r(i,j,4)-r(i,j,3));
i=linspace(10,11,2);
E(i,j,h)=(Q3+2/3*F)*(c(i+1,j,h)-c(i,j,h))+beta*Q3*(c(i,j-1,h)-c(i,j,h))-M*(r(i,j,4)-r(i,j,3));
i=12;
E(i,j,h)=(Q3+2/3*F)*(c(i,j-1,h)-c(i,j,h))+10*beta*Q3*(c(i+1,j,h)-c(i,j,h))-M*(r(i,j,4)-r(i,j,3));
i=13;
E(i,j,h)=Q4*(c(i,j-1,h)-c(i,j,h))+10*beta*Q3*(c(i-1,j,h)-c(i,j,h))-M*(r(i,j,4)-r(i,j,3));
i=linspace(14,15,2);
E(i,j,h)=Q4*(c(i-1,j,h)-c(i,j,h))+beta*Q4*(c(i,j-1,h)-c(i,j,h))-M*(r(i,j,4)-r(i,j,3));
i=16;
E(i,j,h)=F*0+Q4*c(i-1,j,h)-(F+Q4)*c(i,j,h)-M*(r(i,j,4)-r(i,j,3));
%------------------------------------------------------------------
h=5;
j=1;
i=1;
E(i,j,h)=Q1*(c(i,j+1,h)-c(i,j,h))+M*r(i,j,4);
i=linspace(2,3,2);
E(i,j,h)=Q1*(c(i-1,j,h)-c(i,j,h))+beta*Q1*(c(i,j+1,h)-c(i,j,h))+M*r(i,j,4);
i=4;
E(i,j,h)=Q1*c(i-1,j,h)+F*c(i+1,j,h)-(Q1+F)*c(i,j,h)+10*beta*Q1*(c(i+1,j,h)-c(i,j,h))+M*r(i,j,4);
i=5;
E(i,j,h)=(Q2+2/3*F)*(c(i+1,j,h)-c(i,j,h))+10*beta*Q1*(c(i-1,j,h)-c(i,j,h))+M*r(i,j,4);
i=linspace(6,7,2);
E(i,j,h)=(Q2+2/3*F)*(c(i+1,j,h)-c(i,j,h))+beta*Q2*(c(i,j+1,h)-c(i,j,h))+M*r(i,j,4);
i=8;
E(i,j,h)=(Q2+2/3*F)*(c(i,j+1,h)-c(i,j,h))+beta*Q2*(c(i+1,j,h)-c(i,j,h))+M*r(i,j,4);
%第二个循环与第三个循环之间物料交换采用不同的流量以上一个循环为标准
i=9;
E(i,j,h)=Q3*(c(i,j+1,h)-c(i,j,h))+beta*Q2*(c(i-1,j,h)-c(i,j,h))+M*r(i,j,4);
i=linspace(10,11,2);
E(i,j,h)=Q3*(c(i-1,j,h)-c(i,j,h))+beta*Q3*(c(i,j+1,h)-c(i,j,h))+M*r(i,j,4);
i=12;
E(i,j,h)=Q3*c(i-1,j,h)+F*c(i+1,j,h)-(Q3+F)*c(i,j,h)+10*beta*Q3*(c(i+1,j,h)-c(i,j,h))+M*r(i,j,4);
i=13;
E(i,j,h)=(Q4+2/3*F)*(c(i+1,j,h)-c(i,j,h))+10*beta*Q3*(c(i-1,j,h)-c(i,j,h))+M*r(i,j,4);
i=linspace(14,15,2);
E(i,j,h)=(Q4+2/3*F)*(c(i+1,j,h)-c(i,j,h))+beta*Q4*(c(i,j+1,h)-c(i,j,h))+M*r(i,j,4);
i=16;
E(i,j,h)=(Q4+2/3*F)*(c(i,j+1,h)-c(i,j,h))+M*r(i,j,4);
%------------------------------------------------------------------------
j=2;
i=1;
E(i,j,h)=Q1*(c(i,j+1,h)-c(i,j,h))+beta*Q1*(c(i+1,j,h)-c(i,j,h))+M*r(i,j,4);
i=2;
E(i,j,h)=Q1*(c(i,j+1,h)-c(i,j,h))+beta*Q1*(c(i-1,j,h)+c(i,j-1,h)-2*c(i,j,h))+M*r(i,j,4);
i=3;
E(i,j,h)=Q1*c(i-1,j,h)+1/3*F*c(i+1,j,h)-(Q1+1/3*F)*c(i,j,h)+beta*Q1*(c(i,j-1,h)+c(i+1,j,h)-2*c(i,j,h))+M*r(i,j,4);
i=4;
E(i,j,h)=(Q1+F)*(c(i,j-1,h)-c(i,j,h))+10*beta*Q1*(c(i+1,j,h)-c(i,j,h))+beta*Q1*(c(i-1,j,h)-c(i,j,h))+M*r(i,j,4);
i=5;
E(i,j,h)=(Q2-1/3*F)*c(i,j-1,h)+1/3*F*c(i+1,j,h)-Q2*c(i,j,h)+10*beta*Q1*(c(i-1,j,h)-c(i,j,h))+beta*Q2*(c(i+1,j,h)-c(i,j,h))+M*r(i,j,4);
i=6;
E(i,j,h)=(Q2+1/3*F)*(c(i+1,j,h)-c(i,j,h))+beta*Q2*(c(i,j-1,h)+c(i-1,j,h)-2*c(i,j,h))+M*r(i,j,4);
i=7;
E(i,j,h)=(Q2+1/3*F)*(c(i,j+1,h)-c(i,j,h))+beta*Q2*(c(i,j-1,h)+c(i+1,j,h)-2*c(i,j,h))+M*r(i,j,4);
i=8;
E(i,j,h)=(Q2+2/3*F)*(c(i,j+1,h)-c(i,j,h))+beta*Q2*(c(i-1,j,h)-c(i,j,h))+beta*Q2*(c(i+1,j,h)-c(i,j,h))+M*r(i,j,4);
i=9;
E(i,j,h)=Q3*(c(i,j+1,h)-c(i,j,h))+beta*Q2*(c(i-1,j,h)-c(i,j,h))+beta*Q3*(c(i+1,j,h)-c(i,j,h))+M*r(i,j,4);
i=10;
E(i,j,h)=Q3*(c(i,j+1,h)-c(i,j,h))+beta*Q3*(c(i-1,j,h)+c(i,j-1,h)-2*c(i,j,h))+M*r(i,j,4);
i=11;
E(i,j,h)=Q3*c(i-1,j,h)+1/3*F*c(i+1,j,h)-(Q3+1/3*F)*c(i,j,h)+beta*Q3*(c(i,j-1,h)+c(i+1,j,h)-2*c(i,j,h))+M*r(i,j,4);
i=12;
E(i,j,h)=(Q3+F)*(c(i,j-1,h)-c(i,j,h))+beta*Q3*(c(i-1,j,h)-c(i,j,h))+10*beta*Q3*(c(i+1,j,h)-c(i,j,h))+M*r(i,j,4);
i=13;
E(i,j,h)=(Q4-1/3*F)*c(i,j-1,h)+1/3*F*c(i+1,j,h)-Q4*c(i,j,1)+10*beta*Q3*(c(i-1,j,h)-c(i,j,h))+beta*Q4*(c(i+1,j,h)-c(i,j,h))+M*r(i,j,4);
i=14;
E(i,j,h)=(Q4+1/3*F)*(c(i+1,j,h)-c(i,j,h))+beta*Q4*(c(i-1,j,h)+c(i,j-1,h)-2*c(i,j,h))+M*r(i,j,4);
i=15;
E(i,j,h)=(Q4+1/3*F)*(c(i,j+1,h)-c(i,j,h))+beta*Q4*(c(i+1,j,h)+c(i,j-1,h)-2*c(i,j,h))+M*r(i,j,4);
i=16;
E(i,j,h)=(Q4+2/3*F)*(c(i,j+1,h)-c(i,j,h))+beta*Q4*(c(i-1,j,h)-c(i,j,h))+M*r(i,j,4);
%-----------------------------------------------------------------------
j=3;
i=1;
E(i,j,h)=1/3*F*c(i+1,j,h)+(Q1-1/3*F)*c(i,j+1,h)-Q1*c(i,j,h)+beta*Q1*(c(i+1,j,h)-c(i,j,h))+M*r(i,j,4);
i=2;
E(i,j,h)=(Q1+1/3*F)*(c(i+1,j,h)-c(i,j,h))+beta*Q1*(c(i,j+1,h)+c(i-1,j,h)-2*c(i,j,h))+M*r(i,j,4);
i=3;
E(i,j,h)=(Q1+1/3*F)*(c(i,j-1,h)-c(i,j,h))+beta*Q1*(c(i+1,j,h)+c(i,j+1,h)-2*c(i,j,h))+M*r(i,j,4);
i=4;
E(i,j,h)=(Q1+2/3*F)*(c(i,j-1,h)-c(i,j,h))+beta*Q1*(c(i-1,j,h)-c(i,j,h))+10*beta*Q1*(c(i+1,j,h)-c(i,j,h))+M*r(i,j,4);
i=5;
E(i,j,h)=Q2*(c(i,j-1,h)-c(i,j,1))+10*beta*Q1*(c(i-1,j,h)-c(i,j,h))+beta*Q2*(c(i+1,j,h)-c(i,j,h))+M*r(i,j,4);
i=6;
E(i,j,h)=Q2*(c(i,j-1,h)-c(i,j,h))+beta*Q2*(c(i-1,j,h)+c(i,j+1,h)-2*c(i,j,h))+M*r(i,j,4);
i=7;
E(i,j,h)=Q2*c(i-1,j,h)+1/3*F*c(i+1,j,h)-(Q2+1/3*F)*c(i,j,h)+beta*Q2*(c(i,j+1,h)+c(i+1,j,h)-2*c(i,j,h))+M*r(i,j,4);
i=8;
E(i,j,h)=(Q2+F)*(c(i,j+1,h)-c(i,j,h))+beta*Q2*(c(i-1,j,h)-c(i,j,h))+beta*Q2*(c(i+1,j,h)-c(i,j,h))+M*r(i,j,4);
i=9;
E(i,j,h)=(Q3-1/3*F)*c(i,j+1,h)+1/3*F*c(i+1,j,h)-Q3*c(i,j,h)+beta*Q2*(c(i-1,j,h)-c(i,j,h))+beta*Q3*(c(i+1,j,h)-c(i,j,h))+M*r(i,j,4);
i=10;
E(i,j,h)=(Q3+1/3*F)*(c(i+1,j,h)-c(i,j,h))+beta*Q3*(c(i,j+1,h)+c(i-1,j,h)-2*c(i,j,h))+M*r(i,j,4);
i=11;
E(i,j,h)=(Q3+1/3*F)*(c(i,j-1,h)-c(i,j,h))+beta*Q3*(c(i,j+1,h)+c(i+1,j,h)-2*c(i,j,h))+M*r(i,j,4);
i=12;
E(i,j,h)=(Q3+2/3*F)*(c(i,j-1,h)-c(i,j,h))+beta*Q3*(c(i-1,j,h)-c(i,j,h))+10*beta*Q3*(c(i+1,j,h)-c(i,j,h))+M*r(i,j,4);
i=13;
E(i,j,h)=Q4*(c(i,j-1,h)-c(i,j,h))+10*beta*Q3*(c(i-1,j,h)-c(i,j,h))+beta*Q4*(c(i+1,j,h)-c(i,j,h))+M*r(i,j,4);
i=14;
E(i,j,h)=Q4*(c(i,j-1,h)-c(i,j,h))+beta*Q4*(c(i-1,j,h)+c(i,j+1,h)-2*c(i,j,h))+M*r(i,j,4);
i=15;
E(i,j,h)=Q4*c(i-1,j,h)+1/3*F*c(i+1,j,h)-(Q4+1/3*F)*c(i,j,h)+beta*Q4*(c(i,j+1,h)+c(i+1,j,h)-2*c(i,j,h))+M*r(i,j,4);
i=16;
E(i,j,h)=(Q4+F)*(c(i,j+1,h)-c(i,j,h))+beta*Q4*(c(i-1,j,h)-c(i,j,h))+M*r(i,j,4);
%--------------------------------------------------------------------------
j=4;
i=1;
E(i,j,h)=(Q1+2/3*F)*(c(i+1,j,h)-c(i,j,h))+M*r(i,j,4);%出口处平衡方程
i=linspace(2,3,2);
E(i,j,h)=(Q1+2/3*F)*(c(i+1,j,h)-c(i,j,h))+beta*Q1*(c(i,j-1,h)-c(i,j,h))+M*r(i,j,4);
i=4;
E(i,j,h)=(Q1+2/3*F)*(c(i,j-1,h)-c(i,j,h))+10*beta*Q1*(c(i+1,j,h)-c(i,j,h))+M*r(i,j,4);
i=5;
E(i,j,h)=Q2*(c(i,j-1,h)-c(i,j,h))+10*beta*Q1*(c(i-1,j,h)-c(i,j,h))+M*r(i,j,4);
i=linspace(6,7,2);
E(i,j,h)=Q2*(c(i-1,j,h)-c(i,j,h))+beta*Q2*(c(i,j-1,h)-c(i,j,h))+M*r(i,j,4);
i=8;
E(i,j,h)=Q2*c(i-1,j,h)+F*c(i+1,j,h)-(Q2+F)*c(i,j,h)+beta*Q2*(c(i+1,j,h)-c(i,j,h))+M*r(i,j,4);
i=9;
E(i,j,h)=(Q3+2/3*F)*(c(i+1,j,h)-c(i,j,h))+beta*Q2*(c(i-1,j,h)-c(i,j,h))+M*r(i,j,4);
i=linspace(10,11,2);
E(i,j,h)=(Q3+2/3*F)*(c(i+1,j,h)-c(i,j,h))+beta*Q3*(c(i,j-1,h)-c(i,j,h))+M*r(i,j,4);
i=12;
E(i,j,h)=(Q3+2/3*F)*(c(i,j-1,h)-c(i,j,h))+10*beta*Q3*(c(i+1,j,h)-c(i,j,h))+M*r(i,j,4);
i=13;
E(i,j,h)=Q4*(c(i,j-1,h)-c(i,j,h))+10*beta*Q3*(c(i-1,j,h)-c(i,j,h))+M*r(i,j,4);
i=linspace(14,15,2);
E(i,j,h)=Q4*(c(i-1,j,h)-c(i,j,h))+beta*Q4*(c(i,j-1,h)-c(i,j,h))+M*r(i,j,4);
i=16;
E(i,j,h)=F*0+Q4*c(i-1,j,h)-(F+Q4)*c(i,j,h)+M*r(i,j,4);
%--------------------------------------------------------
h=6;
j=1;
i=1;
E(i,j,h)=Q1*(c(i,j+1,h)-c(i,j,h))-M*(r(i,j,6)-r(i,j,5));
i=linspace(2,3,2);
E(i,j,h)=Q1*(c(i-1,j,h)-c(i,j,h))+beta*Q1*(c(i,j+1,h)-c(i,j,h))-M*(r(i,j,6)-r(i,j,5));
i=4;
E(i,j,h)=Q1*c(i-1,j,h)+F*c(i+1,j,h)-(Q1+F)*c(i,j,h)+10*beta*Q1*(c(i+1,j,h)-c(i,j,h))-M*(r(i,j,6)-r(i,j,5));
i=5;
E(i,j,h)=(Q2+2/3*F)*(c(i+1,j,h)-c(i,j,h))+10*beta*Q1*(c(i-1,j,h)-c(i,j,h))-M*(r(i,j,6)-r(i,j,5));
i=linspace(6,7,2);
E(i,j,h)=(Q2+2/3*F)*(c(i+1,j,h)-c(i,j,h))+beta*Q2*(c(i,j+1,h)-c(i,j,h))-M*(r(i,j,6)-r(i,j,5));
i=8;
E(i,j,h)=(Q2+2/3*F)*(c(i,j+1,h)-c(i,j,h))+beta*Q2*(c(i+1,j,h)-c(i,j,h))-M*(r(i,j,6)-r(i,j,5));
%第二个循环与第三个循环之间物料交换采用不同的流量以上一个循环为标准
i=9;
E(i,j,h)=Q3*(c(i,j+1,h)-c(i,j,h))+beta*Q2*(c(i-1,j,h)-c(i,j,h))-M*(r(i,j,6)-r(i,j,5));
i=linspace(10,11,2);
E(i,j,h)=Q3*(c(i-1,j,h)-c(i,j,h))+beta*Q3*(c(i,j+1,h)-c(i,j,h))-M*(r(i,j,6)-r(i,j,5));
i=12;
E(i,j,h)=Q3*c(i-1,j,h)+F*c(i+1,j,h)-(Q3+F)*c(i,j,h)+10*beta*Q3*(c(i+1,j,h)-c(i,j,h))-M*(r(i,j,6)-r(i,j,5));
i=13;
E(i,j,h)=(Q4+2/3*F)*(c(i+1,j,h)-c(i,j,h))+10*beta*Q3*(c(i-1,j,h)-c(i,j,h))-M*(r(i,j,6)-r(i,j,5));
i=linspace(14,15,2);
E(i,j,h)=(Q4+2/3*F)*(c(i+1,j,h)-c(i,j,h))+beta*Q4*(c(i,j+1,h)-c(i,j,h))-M*(r(i,j,6)-r(i,j,5));
i=16;
E(i,j,h)=(Q4+2/3*F)*(c(i,j+1,h)-c(i,j,h))-M*(r(i,j,6)-r(i,j,5));
%------------------------------------------------------------------------
j=2;
i=1;
E(i,j,h)=Q1*(c(i,j+1,h)-c(i,j,h))+beta*Q1*(c(i+1,j,h)-c(i,j,h))-M*(r(i,j,6)-r(i,j,5));
i=2;
E(i,j,h)=Q1*(c(i,j+1,h)-c(i,j,h))+beta*Q1*(c(i-1,j,h)+c(i,j-1,h)-2*c(i,j,h))-M*(r(i,j,6)-r(i,j,5));
i=3;
E(i,j,h)=Q1*c(i-1,j,h)+1/3*F*c(i+1,j,h)-(Q1+1/3*F)*c(i,j,h)+beta*Q1*(c(i,j-1,h)+c(i+1,j,h)-2*c(i,j,h))-M*(r(i,j,6)-r(i,j,5));
i=4;
E(i,j,h)=(Q1+F)*(c(i,j-1,h)-c(i,j,h))+10*beta*Q1*(c(i+1,j,h)-c(i,j,h))+beta*Q1*(c(i-1,j,h)-c(i,j,h))-M*(r(i,j,6)-r(i,j,5));
i=5;
E(i,j,h)=(Q2-1/3*F)*c(i,j-1,h)+1/3*F*c(i+1,j,h)-Q2*c(i,j,h)+10*beta*Q1*(c(i-1,j,h)-c(i,j,h))+beta*Q2*(c(i+1,j,h)-c(i,j,h))-M*(r(i,j,6)-r(i,j,5));
i=6;
E(i,j,h)=(Q2+1/3*F)*(c(i+1,j,h)-c(i,j,h))+beta*Q2*(c(i,j-1,h)+c(i-1,j,h)-2*c(i,j,h))-M*(r(i,j,6)-r(i,j,5));
i=7;
E(i,j,h)=(Q2+1/3*F)*(c(i,j+1,h)-c(i,j,h))+beta*Q2*(c(i,j-1,h)+c(i+1,j,h)-2*c(i,j,h))-M*(r(i,j,6)-r(i,j,5));
i=8;
E(i,j,h)=(Q2+2/3*F)*(c(i,j+1,h)-c(i,j,h))+beta*Q2*(c(i-1,j,h)-c(i,j,h))+beta*Q2*(c(i+1,j,h)-c(i,j,h))-M*(r(i,j,6)-r(i,j,5));
i=9;
E(i,j,h)=Q3*(c(i,j+1,h)-c(i,j,h))+beta*Q2*(c(i-1,j,h)-c(i,j,h))+beta*Q3*(c(i+1,j,h)-c(i,j,h))-M*(r(i,j,6)-r(i,j,5));
i=10;
E(i,j,h)=Q3*(c(i,j+1,h)-c(i,j,h))+beta*Q3*(c(i-1,j,h)+c(i,j-1,h)-2*c(i,j,h))-M*(r(i,j,6)-r(i,j,5));
i=11;
E(i,j,h)=Q3*c(i-1,j,h)+1/3*F*c(i+1,j,h)-(Q3+1/3*F)*c(i,j,h)+beta*Q3*(c(i,j-1,h)+c(i+1,j,h)-2*c(i,j,h))-M*(r(i,j,6)-r(i,j,5));
i=12;
E(i,j,h)=(Q3+F)*(c(i,j-1,h)-c(i,j,h))+beta*Q3*(c(i-1,j,h)-c(i,j,h))+10*beta*Q3*(c(i+1,j,h)-c(i,j,h))-M*(r(i,j,6)-r(i,j,5));
i=13;
E(i,j,h)=(Q4-1/3*F)*c(i,j-1,h)+1/3*F*c(i+1,j,h)-Q4*c(i,j,1)+10*beta*Q3*(c(i-1,j,h)-c(i,j,h))+beta*Q4*(c(i+1,j,h)-c(i,j,h))-M*(r(i,j,6)-r(i,j,5));
i=14;
E(i,j,h)=(Q4+1/3*F)*(c(i+1,j,h)-c(i,j,h))+beta*Q4*(c(i-1,j,h)+c(i,j-1,h)-2*c(i,j,h))-M*(r(i,j,6)-r(i,j,5));
i=15;
E(i,j,h)=(Q4+1/3*F)*(c(i,j+1,h)-c(i,j,h))+beta*Q4*(c(i+1,j,h)+c(i,j-1,h)-2*c(i,j,h))-M*(r(i,j,6)-r(i,j,5));
i=16;
E(i,j,h)=(Q4+2/3*F)*(c(i,j+1,h)-c(i,j,h))+beta*Q4*(c(i-1,j,h)-c(i,j,h))-M*(r(i,j,6)-r(i,j,5));
%-----------------------------------------------------------------------
j=3;
i=1;
E(i,j,h)=1/3*F*c(i+1,j,h)+(Q1-1/3*F)*c(i,j+1,h)-Q1*c(i,j,h)+beta*Q1*(c(i+1,j,h)-c(i,j,h))-M*(r(i,j,6)-r(i,j,5));
i=2;
E(i,j,h)=(Q1+1/3*F)*(c(i+1,j,h)-c(i,j,h))+beta*Q1*(c(i,j+1,h)+c(i-1,j,h)-2*c(i,j,h))-M*(r(i,j,6)-r(i,j,5));
i=3;
E(i,j,h)=(Q1+1/3*F)*(c(i,j-1,h)-c(i,j,h))+beta*Q1*(c(i+1,j,h)+c(i,j+1,h)-2*c(i,j,h))-M*(r(i,j,6)-r(i,j,5));
i=4;
E(i,j,h)=(Q1+2/3*F)*(c(i,j-1,h)-c(i,j,h))+beta*Q1*(c(i-1,j,h)-c(i,j,h))+10*beta*Q1*(c(i+1,j,h)-c(i,j,h))-M*(r(i,j,6)-r(i,j,5));
i=5;
E(i,j,h)=Q2*(c(i,j-1,h)-c(i,j,1))+10*beta*Q1*(c(i-1,j,h)-c(i,j,h))+beta*Q2*(c(i+1,j,h)-c(i,j,h))-M*(r(i,j,6)-r(i,j,5));
i=6;
E(i,j,h)=Q2*(c(i,j-1,h)-c(i,j,h))+beta*Q2*(c(i-1,j,h)+c(i,j+1,h)-2*c(i,j,h))-M*(r(i,j,6)-r(i,j,5));
i=7;
E(i,j,h)=Q2*c(i-1,j,h)+1/3*F*c(i+1,j,h)-(Q2+1/3*F)*c(i,j,h)+beta*Q2*(c(i,j+1,h)+c(i+1,j,h)-2*c(i,j,h))-M*(r(i,j,6)-r(i,j,5));
i=8;
E(i,j,h)=(Q2+F)*(c(i,j+1,h)-c(i,j,h))+beta*Q2*(c(i-1,j,h)-c(i,j,h))+beta*Q2*(c(i+1,j,h)-c(i,j,h))-M*(r(i,j,6)-r(i,j,5));
i=9;
E(i,j,h)=(Q3-1/3*F)*c(i,j+1,h)+1/3*F*c(i+1,j,h)-Q3*c(i,j,h)+beta*Q2*(c(i-1,j,h)-c(i,j,h))+beta*Q3*(c(i+1,j,h)-c(i,j,h))-M*(r(i,j,6)-r(i,j,5));
i=10;
E(i,j,h)=(Q3+1/3*F)*(c(i+1,j,h)-c(i,j,h))+beta*Q3*(c(i,j+1,h)+c(i-1,j,h)-2*c(i,j,h))-M*(r(i,j,6)-r(i,j,5));
i=11;
E(i,j,h)=(Q3+1/3*F)*(c(i,j-1,h)-c(i,j,h))+beta*Q3*(c(i,j+1,h)+c(i+1,j,h)-2*c(i,j,h))-M*(r(i,j,6)-r(i,j,5));
i=12;
E(i,j,h)=(Q3+2/3*F)*(c(i,j-1,h)-c(i,j,h))+beta*Q3*(c(i-1,j,h)-c(i,j,h))+10*beta*Q3*(c(i+1,j,h)-c(i,j,h))-M*(r(i,j,6)-r(i,j,5));
i=13;
E(i,j,h)=Q4*(c(i,j-1,h)-c(i,j,h))+10*beta*Q3*(c(i-1,j,h)-c(i,j,h))+beta*Q4*(c(i+1,j,h)-c(i,j,h))-M*(r(i,j,6)-r(i,j,5));
i=14;
E(i,j,h)=Q4*(c(i,j-1,h)-c(i,j,h))+beta*Q4*(c(i-1,j,h)+c(i,j+1,h)-2*c(i,j,h))-M*(r(i,j,6)-r(i,j,5));
i=15;
E(i,j,h)=Q4*c(i-1,j,h)+1/3*F*c(i+1,j,h)-(Q4+1/3*F)*c(i,j,h)+beta*Q4*(c(i,j+1,h)+c(i+1,j,h)-2*c(i,j,h))-M*(r(i,j,6)-r(i,j,5));
i=16;
E(i,j,h)=(Q4+F)*(c(i,j+1,h)-c(i,j,h))+beta*Q4*(c(i-1,j,h)-c(i,j,h))-M*(r(i,j,6)-r(i,j,5));
%--------------------------------------------------------------------------
j=4;
i=1;
E(i,j,h)=(Q1+2/3*F)*(c(i+1,j,h)-c(i,j,h))-M*(r(i,j,6)-r(i,j,5));%出口处平衡方程
i=linspace(2,3,2);
E(i,j,h)=(Q1+2/3*F)*(c(i+1,j,h)-c(i,j,h))+beta*Q1*(c(i,j-1,h)-c(i,j,h))-M*(r(i,j,6)-r(i,j,5));
i=4;
E(i,j,h)=(Q1+2/3*F)*(c(i,j-1,h)-c(i,j,h))+10*beta*Q1*(c(i+1,j,h)-c(i,j,h))-M*(r(i,j,6)-r(i,j,5));
i=5;
E(i,j,h)=Q2*(c(i,j-1,h)-c(i,j,h))+10*beta*Q1*(c(i-1,j,h)-c(i,j,h))-M*(r(i,j,6)-r(i,j,5));
i=linspace(6,7,2);
E(i,j,h)=Q2*(c(i-1,j,h)-c(i,j,h))+beta*Q2*(c(i,j-1,h)-c(i,j,h))-M*(r(i,j,6)-r(i,j,5));
i=8;
E(i,j,h)=Q2*c(i-1,j,h)+F*c(i+1,j,h)-(Q2+F)*c(i,j,h)+beta*Q2*(c(i+1,j,h)-c(i,j,h))-M*(r(i,j,6)-r(i,j,5));
i=9;
E(i,j,h)=(Q3+2/3*F)*(c(i+1,j,h)-c(i,j,h))+beta*Q2*(c(i-1,j,h)-c(i,j,h))-M*(r(i,j,6)-r(i,j,5));
i=linspace(10,11,2);
E(i,j,h)=(Q3+2/3*F)*(c(i+1,j,h)-c(i,j,h))+beta*Q3*(c(i,j-1,h)-c(i,j,h))-M*(r(i,j,6)-r(i,j,5));
i=12;
E(i,j,h)=(Q3+2/3*F)*(c(i,j-1,h)-c(i,j,h))+10*beta*Q3*(c(i+1,j,h)-c(i,j,h))-M*(r(i,j,6)-r(i,j,5));
i=13;
E(i,j,h)=Q4*(c(i,j-1,h)-c(i,j,h))+10*beta*Q3*(c(i-1,j,h)-c(i,j,h))-M*(r(i,j,6)-r(i,j,5));
i=linspace(14,15,2);
E(i,j,h)=Q4*(c(i-1,j,h)-c(i,j,h))+beta*Q4*(c(i,j-1,h)-c(i,j,h))-M*(r(i,j,6)-r(i,j,5));
i=16;
E(i,j,h)=F*0+Q4*c(i-1,j,h)-(F+Q4)*c(i,j,h)-M*(r(i,j,6)-r(i,j,5));
%--------------------------------------------------------------------------
E=E(:);
[ 本帖最后由 风花雪月 于 2006-11-12 07:50 编辑 ] 小弟在版上谢过了 这类问题一般都都自己编程计算的，不过收敛性一般都不太高，你可以试试同轮算法之类的
页:
[1]