【程序改错】代码，错误信息已贴出
本帖最后由 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
=ode45(@ode1,,);
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, ...
回复 1 # romanticer 的帖子
已算出，但想问一下为什么，在外面算了diff，带入ode就好了，ode里面不然算是么？ 回复 1 # romanticer 的帖子
diff差分应该是A的size-1.
在dy(2)的点除中应该要求前后的size都一样的，
举个例子
x=1:10; length是10
length(diff(x))；length是9
页:
[1]