# How to write ODE for matrices?

This problem must have been encountered by many people but I could not find a neat solution googling.

Many problems in mathematical physics can be written as a matrix ordinary differential equation (ODE), i.e. as $$\frac{d A}{dt} = f(A)$$ with initial condition $$A(0) = A_0$$ where $$A$$ (and $$f(A)$$) is an $$n\times n$$ matrix. For simplicity I consider the case where the function $$f$$ is not explicitly time dependent, as in the above equation. In many circumstances it is easy to produce the function $$f$$ with matrix operations. For example for the Schroedinger equation one has $$f(A) = -i[H,A]$$ (square brackets represent the commutator).

Here is the question: what is a neat way to pass the ODE above to NDSolve ? In my situations the matrix $$A$$ (and $$f(A)$$) are hermitian, so it suffices to send to NDSolve the $$n$$ real diagonal terms and the $$n(n-1)/2$$ complex terms in the upper diagonal. Essentially I am trying to find an elegant way to extract equations out of a matrix for general $$n$$.

• Why couldn't I insert displayed equations with  not to speak of \begin{align}?
– lcv
Nov 20, 2019 at 21:27

Providing an example is always useful. Here is a made up example:

SeedRandom;
ic = RandomReal[1, {3,3}];
H = (#+ConjugateTranspose[#])&@RandomComplex[{-1-I, 1+I},{3,3}];

f[A_, B_] := -I (A.B - B.A)


Then:

sol = NDSolveValue[{A'[t] == f[A[t], H], A == ic}, A, {t, 0, 1}];


Visualization:

ReImPlot[sol[t], {t, 0, 1}] • Oh wow that's simple. Since when is this possible? This is indeed very neat. However, how can I be sure that the integrator only uses $n^2$ real variables as opposed to the $n^2$ complex variables of the matrix?
– lcv
Nov 20, 2019 at 23:55
• Also, how can you make it work if $f(A) = [A \circ C , A]$ where $C$ is a fixed matrix and $A \circ C$ is pointwise (Hadamard) product?
– lcv
Nov 21, 2019 at 0:27
• @lcv I don't understand your distinction between real and complex variables. NDSolve will use a matrix of variables, whose values can be complex. As for Hadamard product, in Mathematica that is just A C, that is A * C or Times[A, C]. Nov 21, 2019 at 0:31
• Thanks for your answer. I assume numerically it makes a difference weather variables are real or complex, but maybe in Mathematica it's impossible to make this distinction. Even if this is the case, if I know that the matrix is hermitian I don't need to compute the evolution of the (say) lower triangular part, as it can be obtained from the upper part. For sure this should lead to some improvement. Regarding the Hadamard product, I know that it's * in Mathematica but simply translating my equation does not work with NDSolveValue
– lcv
Nov 21, 2019 at 0:37
• To be precise: setting f[A_]:= (A*C).A - A. (A*C) does not work
– lcv
Nov 21, 2019 at 0:38