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

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[
Jmat.Erf[Diagonal[Flatten[funcs1]/Sqrt[4/Pi + 2*funcs2]]]],
Thread[Flatten@funcs1 == Flatten@mu0 /. t -> 0]];
equations2 =
Flatten@Join[
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.

• Your code sample doesn't work, please check it. (The jmat and pi are suspicious, but this seems not to be the whole story. ) – xzczd Feb 25 '17 at 8:25
• Sorry about that, I've sorted them out (and edited the code above) and it seems to be running now. Do you have any thoughts/advice on how to a) take into account that funcs2 is symmetric or b) speed the thing up in general? For n=10 it's very fast but unsuprisingly becomes sluggish for higher n. – StressedOutStudent Feb 25 '17 at 12:28
• Are you sure it works? I think one still needs to modify Thread[funcs1 == mu0 /. t -> 0] to Thread[Flatten@funcs1 == Flatten@mu0 /. t -> 0] – xzczd Feb 25 '17 at 12:34
• It seemed to run but yes thank you, you're right. – StressedOutStudent Feb 25 '17 at 13:04

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];
D[funcs1, t] == -funcs1 +
Jmat.Erf[Diagonal[Flatten[funcs1]/Sqrt[4/Pi + 2*funcs2]]]],
Thread[Flatten@funcs1 == Flatten@mu0 /. t -> 0]} //. symRule;

equations2 =
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"}];

• Thank you for your help. Do you know of a way of taking into account/enforcing that funcs2 is symmetric? – StressedOutStudent Feb 25 '17 at 13:35
• @StressedOutStudent Check my edit. – xzczd Feb 25 '17 at 14:12
• Perfect, thank you very much for your help. – StressedOutStudent Feb 25 '17 at 14:34
• By the way, your code above produces the plot much faster than mine in the question - do you know why that is? – StressedOutStudent Feb 25 '17 at 15:32
• @StressedOutStudent I think it's because my sols contains less functions. Notice I've removed those duplicated ones with DeleteCases[Flatten@UpperTriangularize@funcs2, 0]. You can also check Length@Flatten@sols. – xzczd Feb 25 '17 at 16:07