# How does NDSolve treat matrices?

I am having trouble trying to solve a differential equation which couples the elements of a matrix. Let's say I have the following

NDSolve[{D[rho[t],t] == f[rho, t], rho[0]==rho0}, rho, {t, 0, 1}]


with f[rho, t] some function of the matrix rho and time. If this function just involves built-in operations on matrices (e.g. Dot product with another matrix) everything is fine. For example, I can successfully solve the system with

f[rho_,t_]:=Dot[rho,M]Cos[t]


where M is a random matrix of the same dimensions of rho. The problem arises when I try to access the components of the matrix. When I define the following function

f[rho_,t_]:=Table[rho[[n,m]],{n,1,2},{m,1,2}]Cos[t]


NDSolve gives me the error message

Part specification rho[t][[0,1]] is longer than depth of object.

(and the same for the other indices). I assume that the all the built-in functions manipulate lists in a different way. What am I suppose to do?

P.s.: The real problem is much more complicated than this. I would not be able to explicitly write the equation for each matrix element.

• see can-i-solve-system-of-differential-equation-in-a-matrix-form for similar question Commented May 14, 2022 at 20:52
• M and Rho0 are undefined Commented May 14, 2022 at 21:58
• The line  f[rho_,t_]:=Table[rho[[n,m]],{n,1,2},{m,1,2}]Cos[t] is inconsistent. Commented May 14, 2022 at 23:18
• rho0 = RandomReal[1, {2, 2}]; NDSolve[{D[rho[t], t] == rho[t] Cos[t], rho[0] == rho0}, rho, {t, 0, 1}] for the second f[] — no need to index the matrix rho[t] since multiplication threads over lists. Commented May 14, 2022 at 23:46
• Thank you everyone. I know it does not make sense to index the matrix in this example, but I would expect it to work since it should reproduce the same matrix. The thing is that in the real code my function has to access the elements of the matrix and sum them in a particular way, which apparently generates those kind of errors in NDSolve. Commented May 15, 2022 at 7:03

I think that what you want to do with your broken code :

f[rho_,t_]:=Table[rho[[n,m]],{n,1,2},{m,1,2}]Cos[t]

is something similar to this code :

Table[Sum[Indexed[rho, n] Indexed[M, {n, m}]
, {n, 1, 2}], {m, 1,  2}] If[t < .5, 1, 2]


I use If[...] instead of Cos[t] because it makes the control of the plot below easier.

rho0 = {1, 2}
M = {{2, 0}, {0, 3}}
ClearAll[f];
f[rho_, t_] :=
Table[Sum[Indexed[rho, n] Indexed[M, {n, m}], {n, 1, 2}], {m, 1,
2}] If[t < .5, 1, 2] (* Cos[t] *)
sol = NDSolveValue[{D[rho[t], {t}] == f[rho[t], t], rho[0] == rho0},
rho, {t, 0, 2}]

Plot[{Indexed[sol[t], {1}], Indexed[sol[t], {2}], Exp[2 t],
2 Exp[3 t]}, {t, 0, 1}, PlotLegends -> Automatic
, PlotStyle -> {Directive[AbsoluteThickness[3], Black],
Directive[AbsoluteThickness[3], Red], Directive[Green],
Directive[Blue]}
, PlotRange -> {{0, .8}, {0, 12}}]


• I have obtained this result with a lot of trial and errors. I don't have clear explanations for the whole. Commented May 14, 2022 at 23:20

An alternative to Indexed to deal with matrices is to prevent f from evaluating until a matrix value is substituted for rho[t]:

f[r_?MatrixQ, t_] := r Cos[5 t] + r[[1, 2]] - Tr[r] Exp[t];
rho0 = {{1, 2}, {3, 4}};
sol = NDSolve[{D[rho[t], t] == f[rho[t], t], rho[0] == rho0},
rho, {t, 0, 10}];