This DAE appears to be of index greater than 1.
6 ビュー (過去 30 日間)
古いコメントを表示
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
0 件のコメント
採用された回答
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 件のコメント
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.
その他の回答 (0 件)
参考
カテゴリ
Help Center および File Exchange で Platform and License についてさらに検索
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!