# How to get DAE representation from equations?

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

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@
NDSolveProcessEquations[{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}]

• Have you looked at the section Index Reduction for DAEs in this documentation article? Jan 7 '20 at 2:09
• The plots are irrelevant to your question, no? Jan 7 '20 at 4:04
• It seems the function expression residual is the $f$ in the implicit form $f=0$ used by the IDA method in NDSolve, like MATLAB's ode15i (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. Jan 8 '20 at 16:58

To get the first-order state-space equations you can use an internal function that NonlinearStateSpaceModel uses under the hood. NonlinearStateSpaceModel cannot handle this at the moment because it has algebraic equations, but the internal function has no problem with such equations.

{{f, h, e}, x} = ControlDEqnsnonaffinestatespaceForm[
eqns, {PX[t], PY[t], z[t]}, {}, {PX[t], PY[t], z[t]}, t, #[[1 ;; 2]] &,
DescriptorStateSpace -> True];


The above will give equations of the form $$e.x'==f$$.

$$h$$ is the outputs which I specified as $$\{PX[t], PY[t], z[t]\}$$.

Thread[{PX[t], PY[t], z[t]} -> h]


And the equations are

{e.D[x, t], Table["==", 5], f}//Transpose // TableForm
`

• It is nice. The problem is that the DAE is not index 1, index reduction haven't been performed. (I migth not have been precise in the question, but it is a requirement.) Jan 7 '20 at 22:57