# Implement NDSolve when I have a very efficient way to compute the right-hand-sides of my coupled ODEs

I am trying speed up NDSolve on a set of coupled non-linear second-order ODEs.

x[1]''[t] = f[1][x[1],x[2],....]
x[2]''[t] = f[2][x[1],x[2],....]
:
x[n]''[t].....


No problem here--these can be dropped right into NDSolve.

However, suppose I have very efficient way to compute the f[i] terms all-at-once. Can I tell NDSolve to take advantage of that function?

Here is a concrete "toy" example of what I am asking. (that is, I am asking how to solve a problem of similar character, not this exact problem).

Suppose the force on a particle is proportional to the sum of the lengths to the other particles. DistanceMatrix is very efficient:

forces[dofs:{p1:{x1_, y1_}, p2:{x2_, y2_}, p3:{x3_, y3_}}] :=
Total[DistanceMatrix[dofs]] dofs


Here is an example set of degrees of freedom:

vars = Flatten[Table[{x[i][t], y[i][t]}, {i, 3}]]


And this might be my right-hand-side:

 Thread[D[#, {t, 2}] & /@ vars == Table[Inactivate[forces][vars][[i]], {i, 1, 6}]]


I believe that this would calculate force 6 times and not just once.

Is there a way to "thread" the left-hand-side vector of second derivatives over the function force[..]

Thanks, Craig

You can use the vectorized version of NDSolveValue. Since you don't provide initial conditions, I made some up:

SeedRandom[1];
sol = NDSolveValue[
{
X''[t] == forces[X[t]],
X[0] == RandomReal[1, {3, 2}],
X'[0] == RandomReal[1, {3, 2}]
},
X,
{t,0,1}
];


Visualization:

Plot[
Evaluate @ Table[Indexed[sol[t], {i,j}], {i,3}, {j,2}],
{t,0,1}
]


• Excellent. Thank you Carl Woll! I’ll mark this as solved. – Craig Carter Oct 23 '18 at 21:27
• Future Readers: I was attempting to generalize this to an arbitrary number of degrees of freedom, and I realized I did't understand this completely. Clear[forces]; forces[dofs_]:= Total[DistanceMatrix[dofs]] dofs doesn't work, but this does: Clear[forces]; forces[dofs_?MatrixQ]:= Total[DistanceMatrix[dofs]] dofs. – Craig Carter Oct 24 '18 at 15:41