Why does replacing zeros in the mass matrix by small numbers solve: "DAE appears to be of index greater than 1"?

Hello,

I've a set of stiff differential equations that I wanted to solve using ode15s. At small initial concentrations the system worked fine, no problems. However, when I increased the initial values by a factor 10 the error following error popped up:

This DAE appears to be of index greater than 1.

In order to speed up the simulations I used a Mass Matrix M with zeros on the diagonal line except for the last cell, a simple example:

0 0 0 0 0 0
0 0 0 0 0 0
0 0 0 0 0 0
0 0 0 0 0 0
0 0 0 0 1 0
0 0 0 0 0 1

As stated previously, this worked fine for small inital values but not for larger initial values. From some literature and other questions I saw that by replacing the zeros on the Mass Matrix diagonal line by small values, e.g. 1e-12, the system would work. Which it did, so this is the mass matrix that I currently use:

1e-12 0     0     0     0 0
0     1e-12 0     0     0 0
0     0     1e-12 0     0 0
0     0     0     1e-12 0 0
0     0     0     0     1 0 
0     0     0     0     0 1

So the problem is solved, only I don't know why. What's the reason that the DAE index is greater than 1 when using zeros, but when using really small values it is fine?

Does somebody have an explanation for this?

As an example, the simplified version of my system (where parameter values are stored in .par):

options = odeset('Mass',M,'RelTol',1e-6,'AbsTol',1e-6);
[t,c]=ode15s(@Equations,par.tspan,par.Initial,options,par); 
function dcdt = Equations(c,par)
% Create a storage
dcdt = zeros(size(c));
% Determine Rate
Rate = par.k .* c(1,1:end-1); 
% Determine Reaction Profile
    dcdt(1,1:end-1) = - Rate;
    dcdt(2,1:end-1) = + Rate;
% Determine Liquid Profile
dcdt(1,end) = dcdt(1,par.Nr+1) + par.A * (par.c0 - c(1,end) ) + par.B * ( par.C - c(1,end) );
dcdt(2,end) = dcdt(2,par.Nr+1) + par.A * (par.c0 - c(2,end) ) + par.B * ( par.C - c(2,end) );
end

Kind regards,

Danny

[Torsten]

Since we don't know your DAE system, we can't tell.

[Danny Helwegen]

Thanks for the quick reply Torsten, I've added some more information to the question.

[Torsten]

So which of the above expressions is dcdt(1:4) for which you set zeros in the mass matrix ?

[Danny Helwegen]

I also should have update the mass matrix, my apologies, but the zeros in the mass matrix are set for dcdt(i,1:end-1). The 1's in the mass matrix are for the Liquid profiles at the end, dcdt(i,end). With i being 1 and 2.

[Torsten]

"dcdt" must be returned to ode15s as a column vector, not as a matrix.

So how do you reshape your matrix "dcdt" at the end of "System_Equations_DFR" to a vector of dimension 5x1 ?

I just want to see the expressions in the column vector dcdt in position 1-4 that correspond to the zero rows of the mass matrix.

[Danny Helwegen]

at the end of the script I make it a column vector. The reason for this specific notation is because I combine the results at the end to account for diffusion effects inside my system that only apply on concentrations till the end-1 (c(i,1:end-1)). So at the end of the script I have this:

        A1 = c(1,1:end-1);
        A2 = c(2,1:end-1);
    B1 = A1(:) * par.X1;
    B2 = A2(:) * par.X2;
C = [B1; B2; B3; B4];

% Combine Results
dcdt = dcdt(:) + C;

[Torsten]

Sorry, but it's impossible for me to recover what is written in dcdt(1:4) from the information given. But the elements of dcdt that correspond to the zero rows of the mass matrix should define algebraic equations that ode15s tries to nullify. Is this the case here ?

[Danny Helwegen]

Yes, that's the case

[Torsten]

Then you should try "fsolve" before calling "ode15s" to calculate consistent values for the algebraic variables of your DAE system. This usually makes the ODE-solver work.

[Danny Helwegen]

Okay, but then I have 2 questions.

First, I can't enter my parameter structure into the fsolve function? Or is there at trick for that.

Secondly, how can I insert the variables into the ode15s?

[Torsten]

If extracting the algebraic equations from the complete DAE system is too difficult to handle, use "decic" before calling ODE15S:

https://de.mathworks.com/help/matlab/ref/decic.html

[Danny Helwegen]

When I try to use decic I get the following error:

Dot indexing is not supported for variables of this type.

I tried to call it with:

    t0 = 0;
    y0 = par.InitCond;
    yp0 = [zeros(length(par.InitCond),1)];
    [y0,yp0] = decic(@System_Equations,t0,y0,yp0,yp0,[],par);

Is there a way to feed par. to the system?

[Torsten]

You read that you have to add "yp" as input vector to System_Equations ?

You read that you have to define fixed_y0 and fixed_yp0 ?

After doing the necessary changes, passing "par" to System_Equations can be done like

[y0,yp0] = decic(@(t,y,yp)System_Equations(t,y,yp,par),t0,y0,fixed_y0,yp0,fixed_yp0);

Read carefully about the list of input arguments to decic:

https://de.mathworks.com/help/matlab/ref/decic.html

[Danny Helwegen]

Even if I say that no conditions are fixed, fixed_y0 = fixed_yp0 = [ ], I recieve the following error:

Error using decic>sls (line 170)
Index may be greater than one.

[Torsten]

Fix the differential variables to their initial values and give initial guesses for the rest (algebraic variables and derivatives). This should fix your system.

If decic is not able to compute consistent initial values, I could look over your equations if you post them in a comprehensible way.