2
$\begingroup$

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?

Improve this question
$\endgroup$
9
  • 2
    $\begingroup$ 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. $\endgroup$
    – Carl Woll
    Commented Jun 29, 2021 at 18:33
  • $\begingroup$ I really need to use the matrix formulation. How exactly would I do that with the foo function? $\endgroup$
    – sam wolfe
    Commented Jun 29, 2021 at 18:34
  • 3
    $\begingroup$ use x'[t] == foo[x[t]] $\endgroup$
    – Carl Woll
    Commented Jun 29, 2021 at 18:35
  • 1
    $\begingroup$ 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. $\endgroup$
    – Carl Woll
    Commented Jun 29, 2021 at 18:43
  • 1
    $\begingroup$ 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. $\endgroup$
    – bRost03
    Commented Jun 30, 2021 at 19:19

2 Answers 2

Reset to default
1
$\begingroup$

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 *)
Improve this answer
$\endgroup$
3
$\begingroup$

As of version 14.1, the problem can be easily solved with NonThreadable attribute:

SeedRandom[141];
coeff = ConstantArray[1, {3, 3}] - IdentityMatrix[3];
ini = RandomReal[1, {3, 3}];

SetAttributes[x, NonThreadable]

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

Plot[sol[t], {t, 0, 1}]

enter image description here

Improve this answer
$\endgroup$

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.