4
$\begingroup$

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

enter image description here

what am I doing wrong?

Improve this question
$\endgroup$
6
  • $\begingroup$ 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? $\endgroup$
    – Nasser
    Commented Aug 26, 2022 at 12:59
  • $\begingroup$ @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? $\endgroup$
    – pixis
    Commented Aug 26, 2022 at 13:03
  • $\begingroup$ 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. $\endgroup$
    – Nasser
    Commented Aug 26, 2022 at 13:14
  • $\begingroup$ @Nasser Evaluation order isn't that easy for beginners… $\endgroup$
    – xzczd ♦
    Commented Aug 26, 2022 at 13:28
  • 1
    $\begingroup$ 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: mathematica.stackexchange.com/questions/18393/… $\endgroup$
    – Michael E2
    Commented Aug 26, 2022 at 13:40

3 Answers 3

Reset to default
6
$\begingroup$

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}]
Improve this answer
$\endgroup$
2
  • $\begingroup$ thanks for your answer! this is exactly what I was confused about $\endgroup$
    – pixis
    Commented Aug 26, 2022 at 13:37
  • 1
    $\begingroup$ @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 :) . $\endgroup$
    – xzczd ♦
    Commented Aug 26, 2022 at 13:39
5
$\begingroup$

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

enter image description here

Improve this answer
$\endgroup$
1
  • $\begingroup$ 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! $\endgroup$
    – pixis
    Commented Aug 26, 2022 at 13:37
5
$\begingroup$

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}];
Improve this answer
$\endgroup$
2
  • $\begingroup$ 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. $\endgroup$
    – Michael E2
    Commented Aug 26, 2022 at 13:37
  • $\begingroup$ Thanks! this works for, I'm accepting the other answer only because it came first and does the same thing $\endgroup$
    – pixis
    Commented Aug 26, 2022 at 13:38

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

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