24小时热门版块排行榜

返回列表

月下行云

新虫 (小有名气)

应助: 2 (幼儿园)
金币: 2718.8
散金: 25
红花: 2
帖子: 170
在线: 416.7小时
虫号: 1528187
注册: 2011-12-07
专业: 控制理论与方法

[求助] matlab解DAE遇到的问题

各位同事：
我在解DAE时，用示例1是可以调试出结果的，但是实例2却总是
Error using daeic12 (line 77)
This DAE appears to be of index greater than 1.

Error in ode15s (line 311)
[y,yp,f0,dfdy,nFE,nPD,Jfac] = daeic12(odeFcn,odeArgs,t,ICtype,Mt,y,yp0,f0,...

Error in sunli_Nsolve (line 12)
[t,x]=ode15s(@odefun,[0 10],x0,options);%

各位知道为什么吗？
示例1程序
function dx = rigid(t,x,k,l,m)
dx= zeros(4,1);
dx(1)=-0.2*x(1)+k*x(2)*x(3)+l*0.3*x(1)*x(2);
dx(2)=2*x(1)*x(2)-5*x(2)*x(3)-2*m*x(2)*x(2);
dx(3)=x(1)+x(2)+x(3)-1;
dx(4)=-x(4)-1+x(3);%微分方程组
end
M=[1 0 0 0;0 1 0 0;0 0 0 0;0 0 0 0;];%质量矩阵
options=odeset('mass',M);%对以DAE问题，mass属性必须设置
x0=[0.8;0.1;0.1;0.1];%初值
k=1;l=1;m=1;
[t,x]=ode15s(@rigid,[0 20],x0,options,k,l,m);%这里好像不能使用ode45
figure('numbertitle','off','name','DAE demo—by Matlabsky')
plot(t,x)
legend('x1(t)','x2(t)','x3(t)')

实例2程序：
function dx = odefun(t,x)
dx=zeros(6,1);
l=1.33;k0=19212;k2=0.000560;
c0=180;c1=10600;c2=59830;d1=500;d2=3000;tao=17;
rb=96;Dw=820;ut=76;Dsw=20;hec=2000;
dx(1)=1/c1*((hec-d1)*Dw+(l*x(2)-hec)*Dsw+(d1-l*x(2))*x(6)+k0*rb);%pm
dx(2)=1/c2*((hec-d2)*Dw+(l*x(2)-hec)*Dsw+(d2-l*x(2))*x(6)+k0*rb);%hm
dx(3)=k2*(x(6)*x(5)-x(5)*hec);
dx(4)=x(1)-(x(6)-Dsw)^2/(74554*(x(1)-x(4))^1.3);
dx(5)=(l*x(2)*(x(6)-Dsw)+Dsw*hec)/x(6);
dx(6)=ut*4877.6*x(4)*(-854+13.47*x(4)+0.4*x(5)-0.003*x(4)*x(5));%微分方程组
end
M=[1 0 0 0 0 0;0 1 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];%质量矩阵
options=odeset('mass',M);%对以DAE问题，mass属性必须设置
x0=[25;2800;950;20;2500;850]
[t,x]=ode15s(@odefun,[0 10],x0,options);%这里好像不能使用ode45
figure('numbertitle','off','name','DAE demo—by Matlabsky')
plot(t,x)
legend('x1(t)','x2(t)','x3(t)')
非常疑惑，基本一样，就加了几个参数

回复此楼