I am using NDsolve for my DAE system. I have 11 Equation and 11 Unknown. I used this code for 8 equation and 8 unknown system and it works correctly but for this system I received this error: Unable to reduce the index of the system to 0 or 1. My questions is : 1) Why this code cant solve my equations?(Where is my mistake?) 2) How can I know my Dae index with Mathematica functions before solve them? 3) Can I reduce my DAE equations system index with Mathematica's function and change them to simple ODE and then solve them with NDSolve?(If yes how can I do it?) 4) I need too solve nonlinear DAE with 30 equation and 30 unknown, can I do it with this method of coding and use of NDsolve or I should change my way?

Clear[ g, LHS, RHS, eqns,contraints, vars, SliderCrankSol, θ1, θ2, θ3,AA, BB, CC, DD, EE, FF, GG, HH, T1, T12, T23, T4]; 

vars = {θ1, θ2, θ3, AA, BB, CC, DD, EE,FF, GG, HH};

T1[t_] := 1;
T12[t_] := 1;
T23[t_] := 1;
T4[t_] := 1;

x1[t_] := -Sin[θ1[t]];
y1[t_] := -Cos[θ1[t]];

x2[t_] := -2*Sin[θ1[t]] - Cos[θ2[t]];
y2[t_] := -2*Cos[θ1[t]] - Sin[θ2[t]];

x3[t_] := -a - Sin[θ3[t]];
y3[t_] :=  b- Cos[θ3[t]];


eq1 = (AA[t] - CC[t])*Sin[θ1[t]] + (DD[t] - BB[t])*Cos[θ1[t]] == D[x1[t], {t, 2}];

eq2 = (AA[t] - CC[t])*Cos[θ1[t]] - (DD[t] - BB[t])*Sin[θ1[t]] - 9.81 == D[y1[t], {t, 2}];

eq3 = -(DD[t] + BB[t]) + T1[t] - T12[t] == D[θ1[t], {t, 2}];

eq4 = CC[t]*Sin[θ1[t]] - DD[t]*Cos[θ1[t]] + FF[t]*Cos[θ3[t]] + EE[t]*Sin[θ3[t]] == D[x2[t], {t, 2}];

eq5 = CC[t]*Cos[θ1[t]] + DD[t]*Sin[θ1[t]] - FF[t]*Sin[θ3[t]] + EE[t]*Cos[θ3[t]] - g == D[y2[t], {t, 2}];

eq6 = CC[t]*Cos[θ1[t] + θ2[t]] + DD[t]*Sin[θ1[t] + θ2[t]] + FF[t]*Sin[θ2[t] + θ3[t]] - EE[t]*Cos[θ2[t] + θ3[t]] - T12[t] + T23[t] ==D[θ2[t], {t, 2}];

eq7 = (GG[t] - EE[t])*Sin[θ3[t]] + (HH[t] - FF[t])*Cos[θ3[t]] == D[x3[t], {t, 2}];

eq8 = (GG[t] - EE[t])*Cos[θ3[t]] - (HH[t] - FF[t])*Sin[θ3[t]] -g == D[y3[t], {t, 2}];

eq9 = +(HH[t] + FF[t]) + T4[t] + T23[t] == D[θ3[t], {t, 2}];

eqns = {eq1, eq2, eq3, eq4, eq5, eq6, eq7, eq8, eq9};
constraints = ( {
{-Sin[θ3[t]] + Cos[θ2[t]] + Sin[θ1[t]] == a/2},{Cos[θ3[t]] - Sin[θ2[t]] -Cos[θ1[t]] == b/2}} );

paramsFSC = { a -> 2, b -> 0, g -> 9.81};
{sliderCrank} = NDSolve[{eqns,constraints, θ1[0] == 0.3, θ1'[0] == 0} /. paramsFSC,
vars, {t, 0, 40}, Method -> {"IndexReduction" -> Automatic}];
Row[Plot[Evaluate[#[t] /. sliderCrank], {t, 0, 20}, PlotRange -> All, 
PlotLabel -> #, ImageSize -> Small] & /@ vars]
  • $\begingroup$ In eq5 the \[Theta] are missing. With Method -> {"IndexReduction" -> {Automatic, "ConstraintMethod" -> None}} you will find a solution to the differential equations but of course the constraints are not met. If you specify a complete set of initial conditions you can met the constraints (with increasing error). You can check the validity of the solution for a given time by {eqns, constraints} /. SetPrecision[sliderCrank, 10] /. g -> 9.81 /. t -> 4 Here are the details. $\endgroup$ Feb 17, 2018 at 13:26
  • $\begingroup$ missing of θ in eq5 occurs after I past my code and while I was editing it. my constraints are important and I cant ignore them. in my system θ2 and θ3 are function of θ1 that is clear from my constraints, so my equation can solve with one initial condition for θ1. I write this code according to ((slidercrank)) example of Mathematica help. my main problem is why this error(Unable to reduce the index of the system to 0 or 1) occurs and how can I resolve it? I use this method for 7 DAEequations and 8 unknown with 1 constraint and it works correctly. $\endgroup$
    – Farshid
    Feb 17, 2018 at 14:00

1 Answer 1


In your example a unique solution is defined by the initial conditions and the constraints from the set of functions that solve the differential equations. If you have a lot of equations it seems to be difficult for Mathematica to do index reduction and keep the equations for the constraints. A solution to this problem is to use the constraints for index reduction but not during solving the differential equations via Method -> {"IndexReduction" -> {Automatic, "ConstraintMethod" -> None}}. The unique solution is then defined adding more initial conditions. The price is that the precision of the solution with respect to the constraints is lower. So this needs to be checked.

First determine initial conditions that meet the constraints:

Flatten[constraints] /. paramsFSC;
Solve[% /. t -> 0 /. \[Theta]1[0] -> 0.3, {\[Theta]2[0], \[Theta]3[0]}]
D[%%, t] /. t -> 0 /. \[Theta]1[0] -> 0.3 /. \[Theta]1'[0] -> 0 /.%[[2]];
Solve[% /. t -> 0 /. \[Theta]1[0] -> 0.3, {\[Theta]2'[0], \[Theta]3'[0]}]

enter image description here


{sliderCrank} = NDSolve[{eqns, constraints, \[Theta]1[0] == 0.3, \[Theta]1'[0] == 0, \[Theta]2[0] == 0, \[Theta]3[0] == 0.3, \[Theta]2'[0] == 0, \[Theta]3'[0] == 0} /. paramsFSC,
 vars, {t, 0, 40}, Method -> {"IndexReduction" -> {Automatic, "ConstraintMethod" -> None}}];
Row[Plot[Evaluate[#[t] /. sliderCrank], {t, 0, 20}, PlotRange -> All, PlotLabel -> #, ImageSize -> Small] & /@ vars]

enter image description here

Then check the validity of the solution for a given time by

{eqns, constraints} /. SetPrecision[sliderCrank, 10] /. paramsFSC /. t -> 30

enter image description here

So up to 10 significant digits the differential equations and the constraints are met.

  • $\begingroup$ thank you dear Matthias for your answer. I understand that initial condition in my DAE system is very important for solve them. If I have more accurate initial condition I can gain more accurate answer. My equations are Newton-Euler equations that are very important in Dynamics and now I understand that direct numerical solution for this type of equations is very challenging issue that is related too initial conditions. but if it can be done that I could generate 3 differential equations of motion from this 9 equation, I could solve 3 ODE with 2 constraint that it was better than this way. $\endgroup$
    – Farshid
    Feb 17, 2018 at 21:58
  • $\begingroup$ How can i find 3 ODE from this 9 equation that solve them with 2 constraint? i think this is a simple and basic way for solve newton-Euler equations. First we extract ODE's from Our DAE's and then solve them like simple ODE with Constraint. I think if any one suggest a basic solution in this manner, this topic can be very useful and Fundamental for any body that works on dynamics of systems. Please offer a way that first find governing ODE's and then NDsolve ODE's with constraint and finally achieve algebraic unknown-AA,BB,CC,... -. if it is not simple or obvious I bring up a new topic for it. $\endgroup$
    – Farshid
    Feb 17, 2018 at 22:18
  • $\begingroup$ Interesting. I guess \[Theta]2'[0] == 0, \[Theta]3'[0] == 0 is determined from D[constraints, t]? BTW, Method option seems not to be necessary here, without it the result doesn't have significant change, but a somewhat different result will be obtained with Method -> {"IndexReduction" -> {"Pantelides", "ConstraintMethod" -> "Projection"}}, not sure which one is reliable. $\endgroup$
    – xzczd
    Feb 18, 2018 at 3:21
  • $\begingroup$ I can't solve this code with Method -> {"IndexReduction" -> {"Pantelides", "ConstraintMethod" -> "Projection"}}, when i use it an error occurred. $\endgroup$
    – Farshid
    Feb 18, 2018 at 6:04
  • $\begingroup$ Please check my Suggestion: First step: Eliminate 8 unknown(AA, BB, CC, DD, EE,FF, GG, HH) and get 3 coupled ODE that consist of θ1,θ2,θ3 and their derivative. second step: NDSolve 3 coupled ODE with 2 constraint. Last step: Find (AA, BB, CC, DD, EE,FF, GG, HH) that are function of θ1,θ2,θ3 and their derivatives. I think this 3 step method is reliable and accurate but i cant write code for extract coupled ODE's and... . If my opinion is wrong please say it to me, and if it's true please help me for execute this 3 step method. $\endgroup$
    – Farshid
    Feb 18, 2018 at 6:24

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.