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

Consider the following test case:

mMax = 50;
Ffuncs[r_] = Table[F[n, m][r], {n, mMax}, {m, mMax}];
NDSolve[Flatten@{
}, 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?

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

With[{opts = SystemOptions[]},
InternalWithLocalSettings[
SetSystemOptions["NDSolveOptions" -> "DefaultSolveTimeConstraint" -> 10.],
sol = NDSolve[
Flatten[D[Ffuncs[r], r]] == Flatten[Ffuncs[r].Ffuncs[r]]],
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}  *)

• But why does NDSolve even need to solve anything algebraically if the derivatives are already given to it? Nov 18, 2016 at 18:55
• @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. Nov 18, 2016 at 20:24
• Note that this doesn't work in Mathematica 9 (it continues to complain about 1 second being too small time). Dec 8, 2016 at 6:08
• @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? Dec 8, 2016 at 11:21
• Might be your machine is too fast. Try increasing mMax. Dec 8, 2016 at 12:29

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.