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],....]

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:

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


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

enter image description here

  • $\begingroup$ Excellent. Thank you Carl Woll! I’ll mark this as solved. $\endgroup$ Oct 23 '18 at 21:27
  • $\begingroup$ 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. $\endgroup$ Oct 24 '18 at 15:41

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.