NDSolve - problems solving for coupled vector/matrix equations of high dimension

Question

I'm in a bit of a pickle. I'm trying to solve the following

$$ \frac{\mathrm{d} \mu (t) }{\mathrm{d} t} = -\mu (t)+JE $$ $$\frac{\mathrm{d} C (t) }{\mathrm{d} t} = JXC(t)+C(t)XJ^T-2C(t)+2TI$$

Where $ \mu $ and E are vectors of dimension $n\times 1$ and $C$, $J$ and $X$ are matrices of dimension $n\times n$ (and $I$ is the identity). $J$ is just sampled from a Gaussian whereas $E$ and $X$ are complicated functions of $ \mu $ and $C$'s coefficients. My code pulled up the error

NDSolveValue::ntdv: Cannot solve to find an explicit formula for the derivatives. Consider using the option Method->{"EquationSimplification"->"Residual"}.

So I tried the suggested method and it returned errors NDSolveValue::nlnum: and NDSolveValue::icfail. I haven't got too much experience with Mathematica and so have no idea how to proceed. It's important to note, the matrix $C(t)$ is symmetric which would help things but I don't know how to take this into account. Here's my code:

n = 10;
jint = 2;
temp = 0.5;
Jmat = ReplacePart[
RandomVariate[
NormalDistribution[0, jint/Sqrt[n]], {n, n}], {i_, i_} -> 0];
(* The following are the initial conditions *)
mu0 = RandomVariate[NormalDistribution[0, 1], {1, n}];
cov0 = Array[0 &, {n, n}];

funcs1 = Array[mu[#1, #2][t] &, {n, 1}];
funcs2 = Array[cov[#1, #2][t] &, {n, n}];
equations1 = 
Flatten@Join[
Thread[D[funcs1, t] == -funcs1 + 
   Jmat.Erf[Diagonal[Flatten[funcs1]/Sqrt[4/Pi + 2*funcs2]]]], 
Thread[Flatten@funcs1 == Flatten@mu0 /. t -> 0]];
equations2 = 
Flatten@Join[
Thread[D[funcs2, t] == 
  2*temp*IdentityMatrix[n] - 2*funcs2 + 
   Jmat.DiagonalMatrix[
     Sqrt[2/(2 + Pi*Diagonal[funcs2])]*
      Exp[Diagonal[Flatten[funcs1^2]/(4/Pi + 2*funcs2)]]].funcs2 +
    funcs2.DiagonalMatrix[
     Sqrt[2/(2 + Pi*Diagonal[funcs2])]*
      Exp[Diagonal[
        Flatten[funcs1^2]/(4/Pi + 2*funcs2)]]].Transpose[Jmat]], 
Thread[funcs2 == cov0 /. t -> 0]];

sols = NDSolveValue[{equations1, equations2}, {funcs1, funcs2}, {t, 0,
10}, Method -> {"EquationSimplification" -> "Residual"}]

I've set n=10 but really I need to look at much higher values.

Any help would be enormously appreciated.

Answer 1

This question is strongly related to this one:

Why does NDSolve need to solve for the derivatives if the equations are already explicitly solved?

In short, just add Method -> {"EquationSimplification" -> "Solve"} to NDSolve/NDSolveValue.

When Method -> {"EquationSimplification" -> "Residual"} is set, NDSolve/NDSolveValue will use a DAE solver to solve the equation system. As far as I can tell, the DAE solver of Mathematica is still not strong enough for many cases, so my personal suggestion is to stick on the ODE solver as much as possible.

Also, it's possible to make use of the symmetry of funcs2. We just need to remove duplicated elements of it, with the help of UpperTriangularize, etc.:

(* Former part of the code is the same as yours. *)
symRule = DeleteCases[
   Rule @@ (Flatten@UpperTriangularize@# & /@ {funcs2\[Transpose], funcs2}) // Thread, 
   0 -> 0];
equations1 = {Thread[
     D[funcs1, t] == -funcs1 + 
       Jmat.Erf[Diagonal[Flatten[funcs1]/Sqrt[4/Pi + 2*funcs2]]]], 
    Thread[Flatten@funcs1 == Flatten@mu0 /. t -> 0]} //. symRule;

equations2 = 
  DeleteCases[{Thread /@ 
      Thread[UpperTriangularize /@ (D[funcs2, t] == 
          2*temp*IdentityMatrix[n] - 2*funcs2 + 
           Jmat.DiagonalMatrix[
             Sqrt[2/(2 + Pi*Diagonal[funcs2])]*
              Exp[Diagonal[Flatten[funcs1^2]/(4/Pi + 2*funcs2)]]].funcs2 + 
           funcs2.DiagonalMatrix[
             Sqrt[2/(2 + Pi*Diagonal[funcs2])]*
              Exp[Diagonal[Flatten[funcs1^2]/(4/Pi + 2*funcs2)]]].Transpose[Jmat])], 
     Thread /@ Thread[UpperTriangularize /@ (funcs2 == cov0) /. t -> 0]}, True, 
    Infinity] //. symRule;
sols = NDSolveValue[{equations1, 
    equations2}, {funcs1, DeleteCases[Flatten@UpperTriangularize@funcs2, 0]} // 
    Flatten, {t, 0, 10}, Method -> {"EquationSimplification" -> "Solve"}];