5
$\begingroup$

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[0] == IdentityMatrix[n];
icZ = z[0] == ConstantArray[0, {n, n, n, n}];
NDSolve[{eqY, eqZ, icY, icZ}, {y, z}, {t, 0, 2}]
$\endgroup$
  • $\begingroup$ 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] == 0, z[0] == {4, 5}}, {y, z}, {t, 0, 2}] $\endgroup$ – xzczd Sep 24 at 11:45
  • $\begingroup$ @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. $\endgroup$ – yarchik Sep 24 at 12:23
  • $\begingroup$ Just to confirm, z[t] is a scalar, right? $\endgroup$ – xzczd Sep 24 at 12:31
  • $\begingroup$ @xzczd No, z[t] is {n,n,n,n} tensor. $\endgroup$ – yarchik Sep 24 at 12:32
  • $\begingroup$ Then your i.c. for z is wrong. You need e.g. icZ = z[0] == ConstantArray[0, {3, 3, 3, 3}];. $\endgroup$ – xzczd Sep 24 at 12:33
5
$\begingroup$

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[0] == IdentityMatrix[n];
icZ = z[0] == ConstantArray[0, {n, n, n, n}];
NDSolve[{eqY, eqZ, icY, icZ}, {y, z}, {t, 0, 2}]

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

| improve this answer | |
$\endgroup$
  • $\begingroup$ Yes, this goes very much in the direction of my actual problem. Let me test this approach. $\endgroup$ – yarchik Sep 24 at 18:23
  • $\begingroup$ Yes, to hide the complicated right-hand sides behind a black-box function is useful. $\endgroup$ – yarchik Sep 25 at 10:15

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

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