the application decic is different from the predefined syntax

Hi There,

In the help center, the syntax of decic is

[y0,yp0] = decic(eqs,vars,constraintEqs,t0,y0_est,fixedVars,yp0_est,options)

However, in the “analyze and manipulate differential algebraic equations”, the employment of decic is

decic(f,t0,[0.98;-0.21;zeros(3,1)],[],zeros(5,1),[],opt)

which is different from the predefined syntax.

Answer 1

[Stephen23]

You are looking at the wrong function help.

The text in that example states "Then use the MATLAB decic function..." so you should be looking at the MATLAB function help https://www.mathworks.com/help/matlab/ref/decic.html which has syntax:

[y0_new,yp0_new] = decic(odefun,t0,y0,fixed_y0,yp0,fixed_yp0,options)

which exactly matches the syntax used in the example https://www.mathworks.com/help/symbolic/analyze-and-manipulate-differential-algebraic-equations.html :

[y0,yp0] = decic(f,t0,[0.98;-0.21;zeros(3,1)],[],zeros(5,1),[],opt)

This also makes sense conceptually, because the symbolic toolbox function DAEFUNCTION converts from symbolic equations to a MATLAB function handle (similarly to MATLABFUNCTION does), so we already know that the Symbolic Toolbox function DECIC (which lists its first inputs as symbolic equations and variables) is not the correct function.

You can also use WHICH to check this yourself.

[Tony Cheng]

Hi Stephen,

I am using ode15s to resolve an index-1 dae system in the form M(y)*dot_y = f(y).

The consistent initial values at t=t0 are provided from our previous computation, then the equality M(y)*dot_y = f(y) is definitely satisfied at this time instant. However, the command window gives an error “Error using daeic3 Need a better guess y0 for consistent initial conditions”.

Does this mean this equality is not satisfied? If so, do I need to use decic to get consistent initial values?

Thanks in advance!

[Stephen23]

This is what my AI tool states:

The message does not prove your equation M(y0)*ydot0 = f(y0) is false. It means ode15s’s internal initializer (daeic3) could not confirm or refine a consistent pair (y0, ydot0) within its tolerances and linear algebra checks. That can happen even when your equality holds in theory, for example due to small residuals relative to tolerances, wrong Mass/MStateDependence options, singular or ill-conditioned Jacobians, or because you only supplied y0 and asked the solver to guess ydot0.

What ode15s needs

For index-1 DAEs, ode15s requires consistent initial conditions: y0 and ydot0 such that r = M(t0,y0)*ydot0 - f(t0,y0) is approximately zero (within tolerances). If you provide only y0, MATLAB calls daeic3 to find ydot0 (and sometimes adjust y0). If that Newton solve fails, you see “Need a better guess y0 for consistent initial conditions”.

What to do

If you already have a consistent initial derivative, pass it explicitly.Example:opts = odeset('Mass',Mfun, ... % your M(t,y) or M(y)'MassSingular','yes', ... % DAE'MStateDependence','strong', ... % if M depends on y'InitialSlope',ydot0, ... % your consistent derivative'RelTol',1e-6,'AbsTol',1e-8);[t,y] = ode15s(@f,tspan,y0,opts);

Notes

• Use MStateDependence 'none' if M is constant, otherwise usually 'strong'.
• If M is singular (true DAE), keep MassSingular 'yes'.

If you only have y0 (or the solver still complains), compute consistent initials with decic before integrating.Example:fixy = false(size(y0)); % allow all y components to move unless you must fix somefixyd = false(size(y0)); % allow all derivatives to move unless you must fix some

[y0c,ydot0c] = decic(@f,t0,y0,fixy,[],fixyd, ...
    'Mass',Mfun, ...
    'MStateDependence','strong', ...
    'RelTol',1e-8,'AbsTol',1e-10);
opts = odeset('Mass',Mfun,'MassSingular','yes','MStateDependence','strong', ...
    'InitialSlope',ydot0c,'RelTol',1e-6,'AbsTol',1e-8);
[t,y] = ode15s(@f,tspan,y0c,opts);

decic lets you control which y components are allowed to change to enforce constraints, which often avoids the daeic3 failure.

Quick checks that often resolve the issue

• Verify the residual yourself at t0:r = Mfun(t0,y0)*ydot0 - f(t0,y0);norm(r)If norm(r) is not small relative to AbsTol/RelTol scaling, ode15s is right to complain.
• Confirm Mass and MStateDependence match your problem:
• MStateDependence 'none' for constant M
• 'strong' when M depends on y
• If some rows of M are (near) zero, those correspond to algebraic equations. Ensure y0 satisfies those algebraic constraints. If not, allow decic to adjust those y components (set the corresponding fixy entries to false).
• Check scaling and tolerances. Poorly scaled states can make a tiny absolute residual large in relative terms. Consider rescaling states or tightening/loosening tolerances appropriately.
• Provide good Jacobians if you can (via odeset options like Jacobian) or ensure f and Mass are smooth. Discontinuities around t0 can derail the initializer.
• If M is constant and nonsingular, you may set MassSingular 'no' and even precompute Minv to rewrite the system as ydot = Minv*f, but only if that is numerically safe.

Bottom line

• The error means the initializer could not verify consistency, not that your equation is inherently wrong.
• If you know ydot0 that makes M(y0)*ydot0 = f(y0) hold, pass it via InitialSlope.
• Otherwise, yes: use decic to compute consistent (y0, ydot0) first. That is the standard and reliable workflow for index-1 DAEs with ode15s.

[Tony Cheng]

Hi Stephen,

Many thx for u and ur AI tool. Actually, my odeset is

option_dae_1 = odeset( 'Mass' , mass_matrix_handle , ...

'RelTol' , 1 * 10 ^ ( - 5 ) , ...

'AbsTol' , 1 * 10 ^ ( - 7 ) , ...

'MStateDependence' , 'strong' , ...

'MassSingular' , 'yes' , ...

'Stats' , 'on' , ...

'BDF' , 'off' ) ;

As suggested, I use decic to obtain yp0_new. In my issue, all 66 elements of y0 and the first 48 elements of yp0 are obtained from our previous study and we do not want them be changed, and the last 18 elements of yp0 can be adjusted. So I set

fixed_y0 = ones( 66 , 1 ) ;

fixed_yp0 = [ ones( 48 , 1 ) ; zeros( 18 , 1 ) ] ;

In the command window, it gave me an error (I translated into English is

The specified components of y0 and yp0 cannot exceed 66.) . a little bit puzzled...

[Stephen23]

You ran into a decic bookkeeping rule. AI says:

Let n = length(y0). In your case, n = 66. decic treats the entries of y0 and yp0 together as the unknowns it is allowed to adjust in order to drive the residual r = M(t0,y0)*yp0 - f(t0,y0) to zero. You tell decic which entries are fixed (cannot move) with fixed_y0 and fixed_yp0. For decic to have enough freedom to solve r = 0 (which has n scalar equations), the total number of free entries across y0 and yp0 must be at least n. Equivalently:

sum(fixed_y0) + sum(fixed_yp0) <= n

Your choices fix all y0 (66 fixed) and fix the first 48 entries of yp0 (48 fixed), so

66 + 48 = 114 > 66

which violates that rule. That’s why you see the error about “specified components … cannot exceed 66”.

What you can do

Option A: If your previous study already produced a consistent pair (y0, yp0), skip decic and pass yp0 directly to ode15s:

opts = odeset('Mass' , mass_matrix_handle , ...

'RelTol' , 1e-5 , ...

'AbsTol' , 1e-7 , ...

'MStateDependence', 'strong' , ...

'MassSingular' , 'yes' , ...

'Stats' , 'on' , ...

'BDF' , 'off' , ...

'InitialSlope' , yp0); % your full-length initial derivative

[t,y] = ode15s(@f, tspan, y0, opts);

If that runs, you do not need decic.

Option B: Use decic, but free enough entries so that the total number of free components across y0 and yp0 is at least 66. Practically:

1. Free the algebraic components of y0. If some rows of M are zero (or nearly so), those y components are algebraic and must be allowed to move to satisfy constraints.
2. Free more entries of yp0 (typically the differential components) until the count of free variables reaches at least n.

Example pattern (you need to define which indices are algebraic vs differential for your model):

% Example masks

is_alg = false(66,1); % set true at algebraic indices

is_diff = ~is_alg;

fixed_y0 = is_diff; % fix differential states, free algebraic states

fixed_yp0 = is_alg; % fix algebraic derivatives, free differential derivatives

% Ensure enough freedom:

assert(sum(fixed_y0) + sum(fixed_yp0) <= 66, 'Too many fixed components')

[y0c, yp0c] = decic(@f, t0, y0, fixed_y0, yp0, fixed_yp0, ...

'Mass', mass_matrix_handle, ...

'MStateDependence', 'strong', ...

'RelTol', 1e-8, 'AbsTol', 1e-10);

opts = odeset('Mass', mass_matrix_handle, 'MassSingular','yes', ...

'MStateDependence','strong', 'InitialSlope', yp0c, ...

'RelTol',1e-5, 'AbsTol',1e-7, 'Stats','on', 'BDF','off');

[t,y] = ode15s(@f, tspan, y0c, opts);

Quick self-checks

1. Count of fixed entries:sumfix = sum(fixed_y0) + sum(fixed_yp0); % must be <= 66
2. If you believe your previous (y0, yp0) are consistent, verify:r = mass_matrix_handle(t0, y0) * yp0 - f(t0, y0);norm(r)

If norm(r) is small relative to your tolerances and scaling, then prefer Option A and pass InitialSlope. If ode15s still complains in that case, double-check that Mass, MassSingular, and MStateDependence are set correctly and that dimensions match.

[Torsten]

Given a component i, you can either specify y0(i) or yp0(i) for "decic", but not both - even if y0(i) and yp0(i) are both correct ("consistent"). Thus the sum of prescribed components of [y0;yp0] is not allowed to exceed 66 in your case.

I encourage you to play with the example below:

2*y1' - y2 = 0

y1 + y2 = 0

odefun = @(t,y,yp) [2*yp(1)-y(2); y(1)+y(2)];
t0 = 0;
y0 = [1 -1];
yp0 = [-0.5 0.5];
[y0,yp0] = decic(odefun,t0,y0,[1 0],yp0,[])      % specifying y1(0) only works

y0 = 2×1

1 -1

<mw-icon class=""></mw-icon>

<mw-icon class=""></mw-icon>

yp0 = 2×1

-0.5000 0.5000

<mw-icon class=""></mw-icon>

<mw-icon class=""></mw-icon>

[y0,yp0] = decic(odefun,t0,y0,[0 1],yp0,[])      % specifying y2(0) only works

y0 = 2×1

1 -1

<mw-icon class=""></mw-icon>

<mw-icon class=""></mw-icon>

yp0 = 2×1

-0.5000 0.5000

<mw-icon class=""></mw-icon>

<mw-icon class=""></mw-icon>

[y0,yp0] = decic(odefun,t0,y0,[1 0],yp0,[0 1])   % specifying y1(0) and y2'(0) works  

y0 = 2×1

1 -1

<mw-icon class=""></mw-icon>

<mw-icon class=""></mw-icon>

yp0 = 2×1

-0.5000 0.5000

<mw-icon class=""></mw-icon>

<mw-icon class=""></mw-icon>

[y0,yp0] = decic(odefun,t0,y0,[1 1],yp0,[1 0])   % specifying y1(0), y1'(0) and y2(0) does not work (similar to your setting)

Error using decic (line 44)
Cannot specify more than 2 components of y0 and yp0.

[y0,yp0] = decic(odefun,t0,y0,[1 0],yp0,[1 0])   % specifying y1(0) and y1'(0) does not work 
[y0,yp0] = decic(odefun,t0,y0,[0 1],yp0,[1 0])   % specifying y2(0) and y1'(0) does not work because y1 is a differential variable

[Tony Cheng]

Dear Stephen,

Thanks very much for your detailed explanation. In fact, all the elements of y at t0 have been computed from our previous work, and the first 48 elements of dot_y at t0 are also available. The equality M(y)*dot_y = f(y) at t0 is met.

When I use decic with all the 66 elements of y fixed while no elements of dot_y fixed to generate y0_new and yp0_new,the command window gave me an error saying that please try to release 18 fixed components.

By contrast, when I do not specify yp0_new in odeset, the command window did not show that error, and the numerical accuracy in terms of displacement, velocity, and acceleration seems a little bit ok for me.

So is decic necessary to be used to generate y0_new and yp0_new before ode15s is used?

Thanks in advance!

[Torsten]

It is absolutly necessary to prescribe the differential variables at t = t0 of your DAE system because differential variables need initial conditions. The algebraic variables or time derivatives of variables are a "nice-to-have", but if you don't specify them and the solver manages to start the integration, all should be fine.