6

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}]
Michael E2
  • 235,386
  • 17
  • 334
  • 747
yarchik
  • 18,202
  • 2
  • 28
  • 66
  • 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}] – xzczd 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. – yarchik Sep 24 '20 at 12:23
  • Just to confirm, z[t] is a scalar, right? – xzczd Sep 24 '20 at 12:31
  • @xzczd No, z[t] is {n,n,n,n} tensor. – yarchik Sep 24 '20 at 12:32
  • Then your i.c. for z is wrong. You need e.g. icZ = z[0] == ConstantArray[0, {3, 3, 3, 3}];. – xzczd Sep 24 '20 at 12:33
  • @xzczd Yes, you are right. Now I have to make my minimal example more complicated to illustrate my point. The equation for y[t] matrix has some complicated contraction in it. NDSolve can solve it. But the problem is with my abstract formulation: the TensorContract/TensorProduct is extremely memory hungry $\mathcal{O}(n^8)$. – yarchik Sep 24 '20 at 12:58
  • See also this question and my own DIY attempt to deal with this issue (not totally successful!) – Chris K Sep 24 '20 at 13:03
  • I think in this case Dot, or Compiled combination of Table and Sum will be a better choice. – xzczd Sep 24 '20 at 13:07
  • @xzczd Yes, I agree with you, and implementing this brought me to the original problem. – yarchik Sep 24 '20 at 13:17
  • @ChrisK Thank you, Chris! I will study your post. – yarchik Sep 24 '20 at 13:19
  • @yarchik not sure it covers vector equations – Chris K Sep 24 '20 at 13:43
  • @yarchik Your example been solved in v.12.0.0.0. But I have different examples when NDSolve gives the message Cannot solve to find an explicit formula for the derivatives. Usually I am using Flatten form of equations and CoefficientArrays to get matrix and put equations in an explicit form. – Alex Trounev Sep 24 '20 at 16:25
  • @AlexTrounev Right, this is also my strategy. But to the best of my understanding, as soon as we convert ODEs into an explicit component-wise form, matrix multiplications in the rhs are being replaced by much slower Sums. It is correct? – yarchik Sep 24 '20 at 17:37
  • @yarchik It depends on a problem we solve and version we use. In you case and in versions 12, 12.1 it is working as it is. In my case with variabes of form y[i,j,k,l][t] I need to make additional step with Flatten and CoefficientArrays to procide with NDSolve. – Alex Trounev Sep 25 '20 at 11:26

1 Answers1

5

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.]

Michael E2
  • 235,386
  • 17
  • 334
  • 747