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

Question

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?

Michael E2 · Accepted Answer · 2016-12-08T12:44:03.850

7

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}  *)

edited Dec 08 '16 at 12:44

answered Nov 17 '16 at 22:02

Michael E2

235,386
17
334
747

3

But why does NDSolve even need to solve anything algebraically if the derivatives are already given to it? – Ruslan Nov 18 '16 at 18:55
1

@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. – Michael E2 Nov 18 '16 at 20:24
Note that this doesn't work in Mathematica 9 (it continues to complain about 1 second being too small time). – Ruslan Dec 08 '16 at 06:08
@Ruslan Hmm, the code works for me without the time constraint option, which doesn't seem to have existed in V9: http://i.stack.imgur.com/4gvHh.png -- But maybe it's solving fast enough on my machine? – Michael E2 Dec 08 '16 at 11:21
Might be your machine is too fast. Try increasing mMax. – Ruslan Dec 08 '16 at 12:29
1

@Ruslan Yes, that was it. The workaround Method -> {"EquationSimplification" -> "Solve"} still works. – Michael E2 Dec 08 '16 at 12:43
@MichaelE2 By replacing sol=NDSolve[...] part in your answer in V10, it gives >"The value of SystemOption !("PreemptiveCheckUseThreads") cannot
be modified". – user55777 Mar 05 '19 at 03:44
@user55777 I mentioned SystemOptions::noset below the code already. Is it a problem for you? – Michael E2 Mar 05 '19 at 13:16

score 3 · Answer 2 · edited Apr 13 '17 at 12:55

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.

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

2 Answers2

Linked