Update2: Thank you for all of your answers. I would like to clarify the question as I did a bit of research on DAE representations and their applications in software. There are 2 representations used in MATLAB:
- ode15s uses the semi-explicit form. We need a function which returns the state derivatives (dxdt) and the constraint evaluation (have to be 0). In the M mass matrix, we have to define which equatuions are the constraint equations (left-hand side is 0) and which are the state derivatives. M is an input of the solver, along with the DAE function.
- ode15i uses the implicit form. In this case, the DAE function by itself describes the DAE completily because the state derivatives are not explicit, so we don't need M matrix.
The DAEs must be index 1 or 0 in both cases.
Both forms can be created using MATLAB Symbolic Toolbox. However, the consistent initial conditions function works only with the implicit form.
So my question is: how to get these DAE representations in Mathematica? I prefer to use Mathematica in symbolic computations, but in this case I may have to use MATLAB Symbolic.
Source of pictures: https://en.wikipedia.org/wiki/Differential-algebraic_system_of_equations
Other links: https://www.mathworks.com/help/matlab/math/solve-differential-algebraic-equations-daes.html https://www.mathworks.com/help/symbolic/solve-differential-algebraic-equations.html
Update1: In the documentation I found methods to get the processed variables and equations. However, I still don't know how to use them to get the form that can be solved in MATLAB. The names of the variables are confusing and I'm not sure why there are more variables than equations.
I mentioned NonlinearStateSpaceModel because the output of this function can be used directly in MATLAB, but it is only for ODEs.
ClearAll["Global`*"] px := PX[t] vx := D[PX[t], t] ax := D[D[PX[t], t], t] py := PY[t] vy := D[PY[t], t] ay := D[D[PY[t], t], t] m := 1 g := -10 r := 1 Fx := z[t]*px Fy := z[t]*py + m*g vars := {PX, PY, z} eqns = {r^2 == px^2 + py^2, m*ax == Fx, m*ay == Fy} state = First@ NDSolve`ProcessEquations[{eqns, PX[0] == 1, PY[0] == 0}, vars, t, Method -> {"IndexReduction" -> Automatic}] residual = state["NumericalFunction"]["FunctionExpression"] MatrixForm[residual[[1]]] MatrixForm[residual[[2]]] Original question: I have the following code which describes a simple pendulum using Newton's law of motion. This is a DAE (differential-algebraic equation) with higher order than one. NDSolve[] can solve the DAE with the settings.
Instead of numerical solving, I would like to get an order-reduced symbolic representation of the DAE: something similar like on this page: https://www.mathworks.com/help/matlab/math/solve-differential-algebraic-equations-daes.html
How to do it in Mathematica? NonlinearStateSpaceModel works really well with ODEs, but I cannot find such a nice function for DAEs.
ClearAll["Global`*"] px := PX[t] vx := D[PX[t], t] ax := D[D[PX[t], t], t] py := PY[t] vy := D[PY[t], t] ay := D[D[PY[t], t], t] m := 1 g := -10 r := 1 Fx := z[t]*px Fy := z[t]*py + m*g eqns = {r^2 == px^2 + py^2, m*ax == Fx, m*ay == Fy} sol = NDSolve[{eqns, PX[0] == 1, PY[0] == 0, PX'[0] == 0, PY'[0] == 0}, {PX, PY, z}, {t, 0, 5}, Method -> {"IndexReduction" -> Automatic}] ParametricPlot[Evaluate[{PX[t], PY[t]} /. sol], {t, 0, 5}] ParametricPlot[Evaluate[{PX[t], z[t]} /. sol], {t, 0, 5}] 
residualis the $f$ in the implicit form $f=0$ used by the IDA method inNDSolve, like MATLAB'sode15i(which may or may not use IDA under the hood). I still can't answer the question about why there are more variables than equations. $\endgroup$