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

Question

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

Answer 1

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}
]