1
$\begingroup$

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.

Improve this question
$\endgroup$
4
  • $\begingroup$ Your code sample doesn't work, please check it. (The jmat and pi are suspicious, but this seems not to be the whole story. ) $\endgroup$
    – xzczd
    Feb 25, 2017 at 8:25
  • $\begingroup$ 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. $\endgroup$
    – StressedOutStudent
    Feb 25, 2017 at 12:28
  • $\begingroup$ 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] $\endgroup$
    – xzczd
    Feb 25, 2017 at 12:34
  • $\begingroup$ It seemed to run but yes thank you, you're right. $\endgroup$
    – StressedOutStudent
    Feb 25, 2017 at 13:04

1 Answer 1

Reset to default
0
$\begingroup$

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"}];
Improve this answer
$\endgroup$
7
  • $\begingroup$ Thank you for your help. Do you know of a way of taking into account/enforcing that funcs2 is symmetric? $\endgroup$
    – StressedOutStudent
    Feb 25, 2017 at 13:35
  • $\begingroup$ @StressedOutStudent Check my edit. $\endgroup$
    – xzczd
    Feb 25, 2017 at 14:12
  • $\begingroup$ Perfect, thank you very much for your help. $\endgroup$
    – StressedOutStudent
    Feb 25, 2017 at 14:34
  • $\begingroup$ By the way, your code above produces the plot much faster than mine in the question - do you know why that is? $\endgroup$
    – StressedOutStudent
    Feb 25, 2017 at 15:32
  • $\begingroup$ @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. $\endgroup$
    – xzczd
    Feb 25, 2017 at 16:07

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge that you have read and understand our privacy policy and code of conduct.

Not the answer you're looking for? Browse other questions tagged or ask your own question.