Using Flatten inside NDSolve

Question

I'm using Flatten to contract a rank 4 tensor and a matrix (rank 2 tensor) inside NDSolve and there seems to be error. Here is my code

g = 1;
tbar = 30;
n = 4;
d0 = IdentityMatrix[n];
d = Table[If[i == j - 1, Sqrt[j], 0], {i, 0, n - 1}, {j, 0, n - 1}];
h = N[Outer[Times, Transpose[d].d, d0] + 
    Outer[Times, d0, Transpose[d].d] + g*
     Outer[Times, d + Transpose[d], 
      d + Transpose[d]]];
c0 = Table[If[i == 1 && j == 1, 1, 0], {i, 0, n - 1}, {j, 0, n - 1}];
sol = NDSolve[{c'[t] == -I*
      Flatten[h, {{1}, {3}, {2, 4}}].Flatten[c[t], {1, 2}], 
    c[0] == c0}, c, {t, 0, tbar}];

I get the error

what am I doing wrong?

Flatten[c[t], {1, 2}] is wrong There is nothing to Flatten here. That is what the error is saying basically. Why do you want to Flatten c[t] for? — Nasser, Aug 26 '22 at 12:59
@Nasser I want to contract the 2nd and 4th indices of the tensor h with the 2 indices of the matrix c. This is a valid operation. I say this because if you did this with h and c0 the output is another matrix. So I'm not sure what you mean when you say this is wrong. can you please explain why? — pixis, Aug 26 '22 at 13:03
can you please explain why? Well, what do you expect the result of Flatten[c[t], {1, 2}] to be? c[t] is not a list. If you just type Flatten[c[t], {1, 2}] on its own, you will see the error. — Nasser, Aug 26 '22 at 13:14
@Nasser my confusion was with evaluation order as is clear now with the other answers. — pixis, Aug 26 '22 at 13:39
Related, although it's primarily about allowing evaluation only when the argument is numeric (as opposed to a tensor), but the basic idea is similar: https://mathematica.stackexchange.com/questions/18393/what-are-the-most-common-pitfalls-awaiting-new-users/26037#26037 — Michael E2, Aug 26 '22 at 13:40

xzczd · Accepted Answer · 2022-08-26T13:31:57.210

6

Another question evaluation order. NDSolve is a function doesn't own attributes like HoldAll, HoldFirst, etc, the equation will evaluate before it's passed into NDSolve: at this point the c[t] isn't a tensor yet! So the warning Flatten::fldep pops up. (For a more detailed explanation, you may want to read this answer. )

One possible solution is, as shown by cvgmt, avoid Flatten of c[t], another possible solution is, define a black-box helper function that evaluates only when its argument is a List:

rhsfunc[c_List] := -I*Flatten[h, {{1}, {3}, {2, 4}}] . Flatten[c, {1, 2}]
sol = NDSolveValue[{c'[t] == rhsfunc[c[t]], c[0] == c0}, c, {t, 0, tbar}]

edited Aug 26 '22 at 13:31

answered Aug 26 '22 at 13:26

xzczd

65,995
9
163
468

thanks for your answer! this is exactly what I was confused about – pixis Aug 26 '22 at 13:37
1

@pixis Glad it helps. Actually you don't need to accept that fast, it's OK to wait for a while (a day or even longer) to see if better answer(s) come up :) . – xzczd Aug 26 '22 at 13:39

cvgmt · Answer 2 · 2022-08-26T13:27:28.040

5

I think you need to Flatten[h, {{1, 3}, {2, 4}}] instead of Flatten[h, {{1}, {3}, {2, 4}}], and Flatten[c0, {1, 2}] instead of Flatten[c[t], {1, 2}].

g = 1;
tbar = 30;
n = 4;
d0 = IdentityMatrix[n];
d = Table[If[i == j - 1, Sqrt[j], 0], {i, 0, n - 1}, {j, 0, n - 1}];
h = N[Outer[Times, Transpose[d] . d, d0] + 
    Outer[Times, d0, Transpose[d] . d] + 
    g*Outer[Times, d + Transpose[d], d + Transpose[d]]];
c0 = Table[If[i == 1 && j == 1, 1, 0], {i, 0, n - 1}, {j, 0, n - 1}];
sol = NDSolve[{c'[t] == -I*Flatten[h, {{1, 3}, {2, 4}}] . c[t], 
   c[0] == Flatten[c0, {1, 2}]}, c, {t, 0, tbar}]
c[2] /. sol[[1]]

edited Aug 26 '22 at 13:27

answered Aug 26 '22 at 13:20

cvgmt

72,231
4
75
133

this is a useful way of doing it, I'm accepting the other answer only because it explains better what i was doing wrong. Thank you! – pixis Aug 26 '22 at 13:37

score 5 · Answer 3 · answered Aug 26 '22 at 13:29

You're trying to flatten c[t] while it's a symbol before its value as a tensor during integration is substituted. If you want to have matrix output value, a standard solution is to write your own flatten function that won't evaluate until it is passed an array.

g = 1;
tbar = 30;
n = 4;
d0 = IdentityMatrix[n];
d = Table[If[i == j - 1, Sqrt[j], 0], {i, 0, n - 1}, {j, 0, n - 1}];
h = N[Outer[Times, Transpose[d] . d, d0] + 
    Outer[Times, d0, Transpose[d] . d] + 
    g*Outer[Times, d + Transpose[d], d + Transpose[d]]];
c0 = Table[If[i == 1 && j == 1, 1, 0], {i, 0, n - 1}, {j, 0, n - 1}];
myFlatten[t_?ArrayQ, spec_] := Flatten[t, spec];
sol = NDSolve[{c'[t] == -I*
      myFlatten[h, {{1}, {3}, {2, 4}}] . myFlatten[c[t], {1, 2}], 
    c[0] == c0}, c, {t, 0, tbar}];

One could also use a generic conditional evaluation function like SetAttributes[when, HoldAll]; when[expr_, test_] /; test := expr and use it thus: when[Flatten[c[t], {1, 2}], ArrayQ[c[t]]] to replace Flatten[c[t], {1, 2}] in the OP. — Michael E2, Aug 26 '22 at 13:37
Thanks! this works for, I'm accepting the other answer only because it came first and does the same thing — pixis, Aug 26 '22 at 13:38

Using Flatten inside NDSolve

3 Answers3