This DAE appears to be of index greater than 1.

5 次查看（过去 30 天）
Dursman Mchabe
Dursman Mchabe 2018-7-26
评论： Dursman Mchabe 2018-8-15
Hi everyone, I am trying to solve 16 DAEs. The code can be seen on the attachment.
if true
% code
function simultaneousEquations
%%EQUATIONS
%dy(3)/dt = 1/A*(B*C-B*y(3))–((y(3)*D*E-F*y(2))/(1/G)+(F/((1+ (H*y(4))/(I*y(5)))*J))
%dy(7)/dt = 1/A*(-B*y(7))–(K*(1+(H*y(4))/(MM*y(8)))(y(7)*D*E/L–y(9)))
%dy(5)/dt = ((y(3)*D*E-F*y(2))/(1/G)+(F/((1+(H*y(4))/(I*y(5)))*J) ) - (0.162*exp(5153/E)*(((y(4)*y(11))/N) - 1)*(O/((y(4)*y(11)) /N)))
%dy(8)/dt = (K*(1+ (H*y(4))/(MM*y(8)))(y(7)* %D*E/L – y(9)))-(-P*Q*R*y(13)*y(14)*(1-(S*y(14))/(1+S*y(14))))
%dy(15) /dt = (-P*Q*R*y(13)*y(14) *(1-(S*y(14))/(1+S*y(14))))- (0*162*exp(-5153/E)*(((y(4)*y(11))/N)-%1*(O/((y(4)*y(11))/N)))
%dy(13)/dt = -y(13)*(-P*Q*R*y(13)*y(14) *(1-(S*y(14))/(1+S*y(14))))*R/T
%dy(16)/dt = (-P*Q*R*y(13)*y(14) *(1- (S*y(14))/(1+S*y(14))))*Z/AA
% y(14) + 2*y(4) - ((y(5)*W*y(14))/(y(14)^2 + W*y(14) + W*X))- %2*((y(5)*W*X)/(y(14)^2 + W*y(14) + W*X)) – ((y(8)*U*y(14))/(y(14)^2 + U*y(14) + U*V)) – 2*((y(8)*U*V-)/(y(14)^2 + U*y(14) + U*V))- Y/y(14) = 0
%U = y(14)*y(6)/y(9)
%V = y(14)*y(10)/y(6)
%W = y(14)*y(1)/y(2)
%X = y(14)*y(11)/y(1)
%Y = y(14)*y(12)
% y(5) = y(2) + y(1) + y(11)
% y(8) = y(9) + y(6) + y(10)
% y(15) = y(9) + y(6) + y(10)
%% INITIAL VALUES
y0 = zeros(16,1); y0(2)= 1.92e-6; y0(3)= 1.7599e-2; y0(4)= 4.879e-3; y0(5)= 1.4e1; y0(7)= 1.336e-4; y0(8)= 4.879e-3; y0(9)= 6.971e-5; y0(11)= 1.238e1; y0(13)= 48.624; y0(14)= 7.413e-6; y0(1)= 1.615; y0(6)= 4.767; y0(10)= 4.212e-5; y0(12)= 1.349e-6; y0(15)= 4.879e-3; y0(16)= 0;
%% PARAMETER VALUES
A = 1.5e-6; B = 1.66667e-5; C = 6.51332e-2; D = 8.314; E = 323.15; F = 149; G = 4.14e-6; H = 1.39e-9; I = 2.89e-9; J = 8.4e-4; K = 9.598e-4; L = 5.15e+3; MM = 3.53e-9; N = 1.07e-7; O = 10; P = 8.825e-3; Q = 12.54; R = 100.0869; S = 0.84; T = 2703; U = 1.7e-3; V =6.55e-8; W = 6.24; X =5.68e-5; Y =5.3e-8; Z = 258.30; AA = 2540;
M = diag([1 1 1 1 1 1 1 0 0 0 0 0 0 0 0 0]); options = odeset('Mass',M,'MassSingular','yes'); tspan = [0 183000]; [t,y] = ode15s(@(ti,yi)revisedModelode(ti,yi,A,B,C,D,E,F,G,H,I,J,K,L,MM,N,O,P,Q,R,S,T,U,V,W,X,Y,Z,AA),tspan,y0,options);
%% FUNCTION
function yp = revisedModelode(t,y,A,B,C,D,E,F,G,H,I,J,K,L,MM,N,O,P,Q,R,S,T,U,V,W,X,Y,Z,AA)
yp=[1/A*(B*C-B*y(3))-((y(3)*D*E-F*y(2))/(1/G)+(F/((1+ (H*y(4))/(I*y(5)))*J)))
1/A*(-B*y(7))-(K*(1+(H*y(4))/(MM*y(8)))*(y(7)*D*E/L-y(9)))
((y(3)*D*E-F*y(2))/(1/G)+(F/((1+(H*y(4))/(I*y(5)))*J) )-(0.162*exp(5153/E)*(((y(4)*y(11))/N) - 1)*(O/((y(4)*y(11)) /N))))
(K*(1+ (H*y(4))/(MM*y(8)))*(y(7)* D*E/L-y(9)))-(-P*Q*R*y(13)*y(14)*(1-(S*y(14))/(1+S*y(14))))
(-P*Q*R*y(13)*y(14) *(1+(S*y(14))/(1-S*y(14))))- (0*162*exp(-5153/E)*(((y(4)*y(11))/N)-1*(O/((y(4)*y(11))/N))))
-y(13)*(-P*Q*R*y(13)*y(14) *(1-(S*y(14))/(1+S*y(14))))*(R/T)
(-P*Q*R*y(13)*y(14) *(1- (S*y(14))/(1+S*y(14))))*(Z/AA)
y(14) + 2*y(4) - ((y(5)*W*y(14))/(y(14)^2 + W*y(14) + W*X))-2*((y(5)*W*X)/(y(14)^2 + W*y(14) + W*X))-((y(8)*U*y(14))/(y(14)^2 + U*y(14) + U*V))-2*((y(8)*U*V)/(y(14)^2 + U*y(14) + U*V))- Y/y(14)
U-(y(14)*y(6)/y(9))
V-(y(14)*y(10)/y(6))
W-(y(14)*y(1)/y(2))
X-(y(14)*y(11)/y(1))
Y-(y(14)*y(12))
y(5) - y(2) - y(1) - y(11)
y(8) - y(9) - y(6) - y(10)
y(15) - y(9) - y(6) - y(10)];
end
I get an error message : "This DAE appears to be of index greater than 1." I tried to follow the link on the previous answers, but it is no longer available.
https://www.mathworks.com/matlabcentral/answers/102944-what-is-the-meaning-of-this-dae-appears-to-be-of-index-greater-than-1-using-ode-solvers-for-solvin
https://www.mathworks.com/help/releases/R2007a/techdoc/index.html?/help/releases/R2007a/techdoc/ref/ode23.html
I have also tried to create a sparse matrix following
https://www.mathworks.com/matlabcentral/answers/108173-error-using-daeic12-this-dae-appears-to-be-of-index-greater-than-1-solution-set-m-sparse-m

请先登录，再回答此问题。

采纳的回答

Torsten
Torsten 2018-8-14
You try to solve
dy(1)/dt = 1/A*(B*C-B*y(3))–((y(3)*D*E-F*y(2))/(1/G)+(F/((1+ (H*y(4))/(I*y(5)))*J))
dy(2)/dt = 1/A*(-B*y(7))–(K*(1+(H*y(4))/(MM*y(8)))(y(7)*D*E/L–y(9)))
...
not
dy(3)/dt = 1/A*(B*C-B*y(3))–((y(3)*D*E-F*y(2))/(1/G)+(F/((1+ (H*y(4))/(I*y(5)))*J))
dy(7)/dt = 1/A*(-B*y(7))–(K*(1+(H*y(4))/(MM*y(8)))(y(7)*D*E/L–y(9)))
...
  3 个评论
显示 1更早的评论
Torsten
Torsten 2018-8-15
编辑：Torsten 2018-8-15
function main
%%EQUATIONS
%d(y(1))/dt = 1/1.5e-6*(1.67e-5*6.51e-2-1.67e-5*(y(1)))(((y(1))*8.314*323.15-149*(y(8)))/(1/4.14e-6)+(149/((1+ (1.39e-9*(y(9)))/(2.89e-9*(y(3))))*8.4e-4))
%d(y(2))/dt = 1/1.5e-6*(-1.67e-5*(y(2)))(9.6e-4*(1+(1.39e-9*(y(9)))/(3.53e-9*(y(4))))((y(2))*8.314*323.15/5.15e3(y(10))))
%d(y(3))/dt = (((y(1))*8.314*323.15-149*(y(8)))/(1/4.14e-6)+(149/((1+(1.39e-9*(y(9)))/(2.89e-9*(y(3))))*8.4e-4) ) - (0.162*exp(5153/323.15)*((((y(9))*(y(11)))/1.1e-7) - 1)*(10/(((y(9))*(y(11))) /1.1e-7)))
%d(y(4))/dt = (9.6e-4*(1+ (1.39e-9*(y(9)))/(3.53e-9*(y(4))))((y(2))* %8.314*323.15/5.15e3 (y(10))))-(-8.825e-3*12.54*100.0869*(y(6))*(y(12))*(1-(0.84*(y(12)))/(1+0.84*(y(12)))))
%d(y(5)) /dt = (-8.825e-3*12.54*100.0869*(y(6))*(y(12))*(1-(0.84*(y(12)))/(1+0.84*(y(12)))))- (0*162*exp(-5153/323.15)*((((y(9))*(y(11)))/1.1e-7)-%1*(10/(((y(9))*(y(11)))/1.1e-7)))
%d(y(6))/dt = -(y(6))*(-8.825e-3*12.54*100.0869*(y(6))*(y(12)) *(1-(0.84*(y(12)))/(1+0.84*(y(12)))))*100.0869/2703
%d(y(7))/dt = (-8.825e-3*12.54*100.0869*(y(6))*(y(12)) *(1- (0.84*(y(12)))/(1+0.84*(y(12)))))*258.30/2540
%(y(12)) + 2*(y(9)) - (((y(3))*6.24*(y(12)))/((y(12))^2 + 6.24*(y(12)) + 6.24*5.68e-5))- 2*(((y(3))*6.24*5.68e-5)/((y(12))^2 + 6.24*(y(12)) + 6.24*5.68e-5)) %(((y(4))*1.7e-3*(y(12)))/((y(12))^2 + 1.7e-3*(y(12)) + 1.7e-3*6.55e-8)) 2*(((y(4))*1.7e-3*6.55e-8-)/((y(12))^2 + 1.7e-3*(y(12)) + 1.7e-3*6.55e-8))- 5.3e-8/(y(12)) = 0
%1.7e-3 = (y(12))*(y(14))/(y(10))
%6.55e-8 = (y(12))*(y(15))/(y(14))
%6.24 = (y(12))*(y(13))/(y(8))
%5.68e-5 = (y(12))*(y(11))/(y(13))
%5.3e-8 = (y(12))*(y(16))
% (y(3)) = (y(8)) + (y(13)) + (y(11))
% (y(4)) = (y(10)) + (y(14)) + (y(15))
% (y(5)) = (y(10)) + (y(14)) + (y(15))
%%INITIAL VALUES
y0 = zeros(16,1);
y0(1)= 0;
y0(2)= 0;
y0(3)= 0;
y0(4)= 0;
y0(5)= 0;
y0(6)= 4.99e-9;
y0(7)= 0;
y0(8)= 0;
y0(9)= 0;
y0(10)= 0;
y0(11)= 0;
y0(12)= 7.413e-6;
y0(13)= 0;
y0(14)= 0;
y0(15)= 0;
y0(16)= 1.349e-6;
M = diag([1 1 1 1 1 1 1 0 0 0 0 0 0 0 0 0]);
options = odeset('Mass',M,'MassSingular','yes');
tspan = [0 183000];
[t,y] = ode15s(@revisedModelode,tspan,y0,options);
end
%%FUNCTION
function yp = revisedModelode(t,y)
yp=zeros(16,1);
yp(1)=(1/1.5e-6)*((1.67e-5*6.51e-2)-((1.67e-5)*(y(1))))-(((y(1))*(8.314)*(323.15)-(149*(y(8))))/((1/4.14e-6)+(149/((1+(1.39e-9*(y(9)))/(2.89e-9*(y(3))))*8.4e-4))));
yp(2)=1/1.5e-6*(-1.67e-5*(y(2)))-(9.6e-4*(1+(1.39e-9*(y(9)))/(3.53e-9*(y(4))))*((y(2))*8.314*323.15/5.15e3-(y(10))));
yp(3)=(((y(1))*8.314*323.15-149*(y(8)))/(1/4.14e-6)+(149/((1+(1.39e-9*(y(9)))/(2.89e-9*(y(3))))*8.4e-4))-(0.162*exp(5153/323.15)*((((y(9))*(y(11)))/1.1e-7)- 1)*(10/(((y(9))*(y(11)))/1.1e-7))));
yp(4)=(9.6e-4*(1+ (1.39e-9*(y(9)))/(3.53e-9*(y(4))))*((y(2))* 8.314*323.15/5.15e3 - (y(10))))-(-8.825e-3*12.54*100.0869*(y(6))*(y(12))*(1-(0.84*(y(12)))/(1+0.84*(y(12)))));
yp(5)=(-8.825e-3*12.54*100.0869*(y(6))*(y(12))*(1-(0.84*(y(12)))/(1+0.84*(y(12)))))- (0*162*exp(-5153/323.15)*((((y(9))*(y(11)))/1.1e-7)-1*(10/(((y(9))*(y(11)))/1.1e-7))));
yp(6)=-(y(6))*(-8.825e-3*12.54*100.0869*(y(6))*(y(12)) *(1-(0.84*(y(12)))/(1+0.84*(y(12)))))*100.0869/2703;
yp(7)=(-8.825e-3*12.54*100.0869*(y(6))*(y(12)) *(1- (0.84*(y(12)))/(1+0.84*(y(12)))))*258.30/2540;
yp(8)=(y(12))+2*(y(9))-(((y(3))*6.24*(y(12)))/((y(12))^2 + 6.24*(y(12))+ 6.24*5.68e-5))-2*(((y(3))*6.24*5.68e-5)/((y(12))^2 + 6.24*(y(12))+6.24*5.68e-5))-(((y(4))*1.7e-3*(y(12)))/((y(12))^2 + 1.7e-3*(y(12)) + 1.7e-3*6.55e-8))-2*(((y(4))*1.7e-3*6.55e-8)/((y(12))^2 + 1.7e-3*(y(12))+ 1.7e-3*6.55e-8))- 5.3e-8/(y(12));
yp(9)=1.7e-3 -(y(12))*(y(14))/(y(10));
yp(10)=6.55e-8 -(y(12))*(y(15))/(y(14));
yp(11)=6.24-(y(12))*(y(13))/(y(8));
yp(12)=5.68e-5-(y(12))*(y(11))/(y(13));
yp(13)=5.3e-8 - (y(12))*(y(16));
yp(14)=(y(3))-(y(8))-(y(13))-(y(11));
yp(15)=(y(4))-(y(10))-(y(14))-(y(15));
yp(16)=(y(5))-(y(10))-(y(14))-(y(15));
end
Now in order to solve your system, the algebraic equations (8)-(16) should uniquely determine y(8)-y(16), given y(1)-y(7). This is not the case. Analyzing your equations, I come to the conclusion that equations (15) and (16) must be used to determine y(14) and y(15). But this is not possible.
Furthermore, subtracting equation (16) from equation (15) leads to y(4)=y(5), but you write different ODEs for y(4) and y(5). This is contraditory.
So summarizing: Recheck your equations for validity.
Best wishes
Torsten.
Dursman Mchabe
Dursman Mchabe 2018-8-15
Thanks a lot Torsten, I will recheck the equations.

请先登录，再进行评论。

更多回答（0 个）

类别

在 Help Center 和 File Exchange 中查找有关 Programming 的更多信息

标签

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by