1

Imagine I want to solve the following matrix system of equations $$ \begin{pmatrix} \dot{x}_{11} & \dot{x}_{12} & \dot{x}_{13}\\ \dot{x}_{21} & \dot{x}_{22} & \dot{x}_{23}\\ \dot{x}_{31} & \dot{x}_{32} & \dot{x}_{33} \end{pmatrix}= \begin{pmatrix} a_{11}x_{11} & a_{12}x_{12} & a_{13}x_{13}\\ a_{21}x_{21} & a_{22}x_{22} & a_{23}x_{23}\\ a_{31}x_{31} & a_{32}x_{32} & a_{33}x_{33} \end{pmatrix} $$ where $x_{ij}\equiv x_{ij}(t)$ and $\{a_{ij}\}$ are some real coefficients. When $\{a_{ij}\}=1$, we can solve this with NDSolveValue as follows

ini = RandomReal[1, {3, 3}];
sol = NDSolveValue[{
   x'[t] == x[t],
   x[0] == ini
   }, x, {t, 0, 1}]

for some random initial conditions ini (also a matrix). Now, how can I include custom coefficients? Specifically, what if the coefficient matrix is such that the diagonal is zero and all the other entries are 1? That is,

$$ \begin{pmatrix} \dot{x}_{11} & \dot{x}_{12} & \dot{x}_{13}\\ \dot{x}_{21} & \dot{x}_{22} & \dot{x}_{23}\\ \dot{x}_{31} & \dot{x}_{32} & \dot{x}_{33} \end{pmatrix}= \begin{pmatrix} 0 & x_{12} & x_{13}\\ x_{21} & 0 & x_{23}\\ x_{31} & x_{32} & 0 \end{pmatrix} $$

I tried setting

coeff = ConstantArray[1, {3, 3}] - IdentityMatrix[3]

followed simply by

ini = RandomReal[1, {3, 3}];
sol = NDSolveValue[{
   x'[t] == coeff * x[t],
   x[0] == ini
   }, x, {t, 0, 1}]

However this doesn't seem to work and I get the error message

enter image description here

which seems to be related with how NDSolveValue is interpreting coeff

enter image description here

I've noticed that, although x[t] within the NDSolve environment can be tretaed as a list or a list of lists, it fails its interpretation in some cases. Any ideas?

sam wolfe
  • 4,663
  • 7
  • 37
  • 2
    why not work with the vector {x11, x12, x13, x21, x22, x23, x31, x32, x33} instead? Then you can use a matrix-vector dot product on the rhs. If you really need to use a matrix formulation, you can define a function foo[x_List] := coeff x and use that on the rhs. – Carl Woll Jun 29 '21 at 18:33
  • I really need to use the matrix formulation. How exactly would I do that with the foo function? – sam wolfe Jun 29 '21 at 18:34
  • 3
    use x'[t] == foo[x[t]] – Carl Woll Jun 29 '21 at 18:35
  • That's great. Any idea why this works in this case? What is x[t] for NDSolve? Is it really a list? I often get confused by this. – sam wolfe Jun 29 '21 at 18:36
  • 1
    Evaluate x'[t] == coeff x[t]. Since x[t] is not a matrix yet, it is treated as a scalar and x[t] is threaded into the matrix. Using foo prevents this early multiplication. – Carl Woll Jun 29 '21 at 18:43
  • I see. I noticed it treats it as a scalar, however, with Dot or . I saw no problem. I've been using it for the case of column systems (the usual one, like you suggested, not matrices). For example, I define a matrix A as the adjacency matrix of a certain graph on which I have my system and simply do A.x[t]. This is still correct, right? – sam wolfe Jun 29 '21 at 19:02
  • @samwolfe yes that's right. You can try it yourself. Compare {{a,b},{c,d}}.x[t] with {{a,b}, {c, d}}*x[t]. – bRost03 Jun 30 '21 at 19:06
  • Interestingly enough, Transpose is also ok (if I want something like $\dot{x}{ij}=x{ji}$). – sam wolfe Jun 30 '21 at 19:11
  • 1
    It is because all of these functions expect their arguments to be Lists. So if you give foo or Dot or Transpose the argument x[t], it won't evaluate and everything will work like you want. – bRost03 Jun 30 '21 at 19:19

1 Answers1

1

As mentioned in the comments, the issue is premature evaluation of the product between coeff and x[t]. Before NDSolve evaluates x[t] with the initial condition, it performs

coeff*x[t] = {{0,x[t],x[t]}, {x[t],0,x[t]},{x[t],x[t],0}}

Then when it substitutes in the initial condition, each of those x[t]s becomes a matrix and you end up with a mess. The idea is to stop the scalar * matrix product happening until x[t] is evaluated to be a matrix. This can be accomplished by using any function which expects a List as an argument and remains unevaluated if the argument is not a list. @CarolWoll suggests in the comments using foo[x_List]:=coeff*x. We can see that

mat={{1,2,3},{4,5,6},{7,8,9}};
coeff*x[t] = {{0,x[t],x[t]},{x[t],0,x[t]},{x[t],x[t],0}}
coeff*x[t] = {{0,x[t],x[t]},{x[t],0,x[t]},{x[t],x[t],0}}/.x[t]->mat =
    {{0,{{1,2,3},{4,5,6},{7,8,9}},{{1,2,3},{4,5,6},{7,8,9}}},{{{1,2,3},
{4,5,6},{7,8,9}},0,{{1,2,3},{4,5,6},{7,8,9}}},{{{1,2,3},{4,5,6},{7,8,9}},
{{1,2,3},{4,5,6},{7,8,9}},0}} (* huge mess *)
foo[x[t]]=foo[x[t]]
foo[mat]={{0,2,3},{4,0,6},{7,8,0}} (* as desired *)
bRost03
  • 2,072
  • 8
  • 14