MATLAB: This DAE appears to be of index greater than 1.

index

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

Best Answer

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