7
$\begingroup$

Consider the following test case:

mMax = 50;
Ffuncs[r_] = Table[F[n, m][r], {n, mMax}, {m, mMax}];
NDSolve[Flatten@{
    Thread[Flatten[D[Ffuncs[r], r]] == Flatten[Ffuncs[r].Ffuncs[r]]],
    Thread[Flatten@Ffuncs[0] == Table[0, {mMax*mMax}]]
    }, Flatten@Ffuncs[r], {r, 0, 2}];

When I set mMax to 49 or less, it works without any problems. But with 50 or higher, I get

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

That was Mathematica 11. Mathematica 9 gives a better message and a better suggestion, which actually works, unlike the suggestion above, which takes forever to complete the solution and eventually crashes the kernel. Here's Mathematica 9's message:

NDSolve::ndsdtc: The time constraint of 1.` seconds was exceeded trying to solve for derivatives, so the system will be treated as a system of differential-algebraic equations. You can use Method->{"EquationSimplification"->"Solve"} to have the system solved as ordinary differential equations. >>

But anyway, I've already gave NDSolve a system of equations explicitly solved for derivatives! Why does it still need to solve for derivatives?

$\endgroup$
5
$\begingroup$

In V10 and higher, you can increase the time limit on Solve with the system option "NDSolveOptions" -> "DefaultSolveTimeConstraint" -> time:

With[{opts = SystemOptions[]},
 Internal`WithLocalSettings[
  SetSystemOptions["NDSolveOptions" -> "DefaultSolveTimeConstraint" -> 10.`],
  sol = NDSolve[
     Flatten@{Thread[
        Flatten[D[Ffuncs[r], r]] == Flatten[Ffuncs[r].Ffuncs[r]]], 
       Thread[Flatten@Ffuncs[0] == Table[0, {mMax*mMax}]]}, 
     Flatten@Ffuncs[r], {r, 0, 2}];,
  SetSystemOptions[opts]
  ]]
(* Spurious message: SystemOptions::noset: ..SystemOption PreemptiveCheckUseThreads..  *)

sol[[1, 1]]
(*  F[1, 1][r] -> InterpolatingFunction[{{0., 2.}}, <>][r]  *)

I'm not sure what takes NDSolve[] so long, but it's much faster than Solve[]:

Solve[Thread[Flatten[D[Ffuncs[r], r]] == Flatten[Ffuncs[r].Ffuncs[r]]], 
   Flatten[D[Ffuncs[r], r]]]; // AbsoluteTiming
(*  {51.1422, Null}  *)
$\endgroup$
  • 2
    $\begingroup$ But why does NDSolve even need to solve anything algebraically if the derivatives are already given to it? $\endgroup$ – Ruslan Nov 18 '16 at 18:55
  • 1
    $\begingroup$ @Ruslan NDSolve constructs a function that represents the "right-hand side" of the system ${\bf x}'(r) = f({\bf x}, r)$. It does that by solving the equations you pass it (I surmise). You don't actually pass derivatives to NDSolve. Each expression has Equal as its head. It has to inspect them somehow to see what's there. I don't understand why it should take so long, though. $\endgroup$ – Michael E2 Nov 18 '16 at 20:24
  • $\begingroup$ Note that this doesn't work in Mathematica 9 (it continues to complain about 1 second being too small time). $\endgroup$ – Ruslan Dec 8 '16 at 6:08
  • $\begingroup$ @Ruslan Hmm, the code works for me without the time constraint option, which doesn't seem to have existed in V9: i.stack.imgur.com/4gvHh.png -- But maybe it's solving fast enough on my machine? $\endgroup$ – Michael E2 Dec 8 '16 at 11:21
  • $\begingroup$ Might be your machine is too fast. Try increasing mMax. $\endgroup$ – Ruslan Dec 8 '16 at 12:29
3
$\begingroup$

It helps to pose this in vector form, see: Vector form using NDSolve

This at least doesn't throw the solving for derivatives error. (tested up to 200..)

mMax = 60;
vf[vals : {_?NumberQ ..}] := 
 Module[{m = Partition[vals, mMax]}, Flatten[m.m]]

vsoln = NDSolveValue[{x'[r] == vf[x[r]], 
     x[0] == Table[RandomReal[.1], {mMax^2}]}, x[r], {r, 0, 1}];

or with the square matix as the unknowm:

mMax = 60;
vf[vals : {{_?NumberQ ..} ..}] := vals.vals
vsoln = NDSolveValue[{x'[r] == vf[x[r]], 
    x[0] == Table[RandomReal[{-1, 1}], {mMax}, {mMax}]}, 
    x[r], {r, 0, 1}];

I threw in a nonzero initial condition so you get a nonzero result. you should verify this is doing the right thing, but I think its correct.

$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

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