# NDSolve: how the initial step of finding explicit derivatives can be avoided?

There are several posts here devoted to resolving the following error message of NDSolve (NDSolve::ntdvdae, NDSolve::ntdv)

Cannot solve to find an explicit formula for the derivatives

The reason for this message is quite obvious: when ordinary differential equations (ODEs) become quite complicated, MA cannot write them in the form

$$\dot{y}_i=f_i(y,t).$$

To the best of my knowledge there are 2 strategies of dealing with it.

1. Formulate ODEs in a vector form (not always possible).
2. Use some Method options (they typically change ODEs to DAEs and employ a simpler solver).

One of the deficiencies of the 2nd route is that it (Implicit Differential-Algebraic IDA) does not work for complex functions. On the other hand, the structure of my ODEs is such that a vector form is not fully possible:

$$\dot{y}(t)=f(y(t),\vec{v}(t),t),\\ \dot{\vec{v}}(t)= M[\vec{v}(t),y(t)].$$

The right hand side (rhs) of the second vector equation can be expressed as some tensor contraction. $$f$$ cannot be expressed as some matrix operation on $$\vec v(t)$$. Even though, the given ODEs have explicit derivatives in the left hand side (lhs), MA tries hard to analyze the rhs, and fails because equations on the right hand side are difficult. So my question is, how can we indicate that the given system of ODEs already has derivatives in an explicit form on the left hand side and urge MA to proceed with a numerical integration.

Just to repeat my question in a simple form. My 1st-order ODEs always have explicit derivatives on the left hand side and a functional of unknown functions and time on the right hand side. How to enforce numerical integration without algebraic preprocessing?

## Edit: an attempt to construct a minimal example

n = 3;
x = RandomReal[1., {n, n, n, n}];
w = RandomReal[1., {n, n}];
w = w + Transpose[w];
Clear[y, z];
eqY = y'[t] == I w.y[t] +
TensorContract[TensorProduct[x, z[t]], {{2, 6}, {3, 7}, {4, 8}}];
eqZ = z'[t] == I z[t] +
Transpose[ y[t].Transpose[y[t].x.y[t], {2, 1, 4, 3}].y[t], {2, 1, 4, 3}];
icY = y == IdentityMatrix[n];
icZ = z == ConstantArray[0, {n, n, n, n}];
NDSolve[{eqY, eqZ, icY, icZ}, {y, z}, {t, 0, 2}]

• Well, without a specific example it's not easy to help. Anyway, "the structure of my ODEs is such that a vector form is not fully possible" this should not be a problem, for example: NDSolve[{y'[t] == 1, z'[t] == {2, 3}, y == 0, z == {4, 5}}, {y, z}, {t, 0, 2}] Sep 24 '20 at 11:45
• @xzczd Thank you, it is not easy to construct a minimal example. But maybe have a look at my edits. ODEs for the matrix function y[t] alone are perfectly solvable. Add an equation for z[t] tensor and everything breaks down. Notice, that I do not get here an error message as in my real example. It can indicate that there is another problem operating here. Sep 24 '20 at 12:23
• Just to confirm, z[t] is a scalar, right? Sep 24 '20 at 12:31
• @xzczd No, z[t] is {n,n,n,n} tensor. Sep 24 '20 at 12:32
• Then your i.c. for z is wrong. You need e.g. icZ = z == ConstantArray[0, {3, 3, 3, 3}];. Sep 24 '20 at 12:33

The current example in the OP does not produce a cannot-solve error message (although it does give a NDSolve::ndsz stiffness/singularity warning). In any case, it sometimes work to hide the complicated right-hand sides behind a black-box function. In this case the NDSolve::ndsz warning goes away.

n = 3;
x = RandomReal[1., {n, n, n, n}];
w = RandomReal[1., {n, n}];
w = w + Transpose[w];
Clear[y, z];
rhsY[x_?ArrayQ, y_?ArrayQ, z_?ArrayQ] :=
I w.y + TensorContract[
TensorProduct[x, z], {{2, 6}, {3, 7}, {4, 8}}];
eqY = y'[t] == rhsY[x, y[t], z[t]];
rhsZ[x_?ArrayQ, y_?ArrayQ, z_?ArrayQ] :=
I z + Transpose[y.Transpose[y.x.y, {2, 1, 4, 3}].y, {2, 1, 4, 3}];
eqZ = z'[t] == rhsZ[x, y[t], z[t]];
icY = y == IdentityMatrix[n];
icZ = z == ConstantArray[0, {n, n, n, n}];
NDSolve[{eqY, eqZ, icY, icZ}, {y, z}, {t, 0, 2}]


[V12.1.1, Mac, in case it matters.]

• Yes, this goes very much in the direction of my actual problem. Let me test this approach. Sep 24 '20 at 18:23
• Yes, to hide the complicated right-hand sides behind a black-box function is useful. Sep 25 '20 at 10:15