|
马上注册,结交更多好友,享用更多功能,让你轻松玩转社区。
您需要 登录 才可以下载或查看,没有账号?我要加入
x
本帖最后由 romanticer 于 2012-3-20 16:35 编辑
试了一下,还是不行啊
function dy=ode1(t,y)
dy=zeros(2,1);
ps=3.03*10^-9; %partial pressure of SO2, atm
Hs=1.24; %Henry constant, M/atm
Ks1=1.29*10^-2; %Ks1
Ks2=6.014*10^-8;%Ks2
po=5*10^-9; %Po3
Ho=9.4*10^-3;
pnh=1.87*10^-9;
Hnh=62;
Knh=1.709*10^-5;
Kw=10^-14;
pho=0.365*10^-9;
Hho=7.1*10^4;
phn=8.54*10^-18;
Hhn=2.1*10^5;
Khn=15.4;
k0=2.4*10^4;
k1=3.7*10^5;
k2=1.5*10^9;
kp1=7.47*10^7;
kp2=13;
ks3=10^3;
ks4=1.02*10^-2;
A=diff(((1+pnh*Hnh*Knh/Kw)*y(2)-(Kw+ps*Hs*Ks1+phn*Hhn*Khn)/y(2)-2*ps*Hs*Ks1*Ks2/(y(2)^2))/(1/(y(2)/ks3+1+ks4/y(2))+2/(y(2)^2/ks3*ks4+y(2)/ks4+1)))
dy(1)=(k0*ps*Hs+k1*ps*Hs*Ks1/y(2)+k2*ps*Hs*Ks1*Ks2/(y(2)^2))*po*Ho+kp1*y(2)*pho*Hho*ps*Hs*Ks1/y(2)/(1+kp2*y(2));
dy(2)=((k0*ps*Hs+k1*ps*Hs*Ks1/y(2)+k2*ps*Hs*Ks1*Ks2/(y(2)^2))*po*Ho+kp1*y(2)*pho*Hho*ps*Hs*Ks1/y(2)/(1+kp2*y(2)))/A;
end
[t,y]=ode45(@ode1,[0 3600],[0 10^-6.17]);
plot(t,y(:,1),t,y(:,2))
A =
[]
??? Error using ==> mldivide
Matrix dimensions must agree.
Error in ==> ode1 at 29
dy(2)=((k0*ps*Hs+k1*ps*Hs*Ks1/y(2)+k2*ps*Hs*Ks1*Ks2/(y(2)^2))*po*Ho+kp1*y(2)*pho*Hho*ps*Hs*Ks1/y(2)/(1+kp2*y(2)))/A;
Error in ==> odearguments at 98
f0 = feval(ode,t0,y0,args{:}); % ODE15I sets args{1} to yp0.
Error in ==> ode45 at 172
[neq, tspan, ntspan, next, t0, tfinal, tdir, y0, f0, odeArgs, odeFcn, ...
|
|