This DAE appears to be of index greater than 1.
6 vues (au cours des 30 derniers jours)
Afficher commentaires plus anciens
Dursman Mchabe
le 26 Juil 2018
Commenté : Dursman Mchabe
le 15 Août 2018
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 commentaires
Réponse acceptée
Torsten
le 14 Août 2018
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 commentaires
Torsten
le 15 Août 2018
Modifié(e) : Torsten
le 15 Août 2018
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.
Plus de réponses (0)
Voir également
Catégories
En savoir plus sur Platform and License dans Help Center et File Exchange
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!