This DAE appears to be of index greater than 1.

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

Answer 1

[Torsten]

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)))
...

[Dursman Mchabe]

Thanks a lot for pointing me in the right direction. I have re-arranged the equations, then solve for: 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))) ... Instead of

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))) ... However the size of y is 16x1 double, whereas for t is 0. I therefore get the error " Error using vertcat Dimensions of arrays being concatenated are not consistent." I don't know how to make the dimensions consistent. The new code is on the link below: https://www.dropbox.com/s/7sz6gli891k0u4l/simultaneousEquations1.m?dl=0

[Torsten]

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]

Thanks a lot Torsten, I will recheck the equations.