2

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

Carl Woll
  • 130,679
  • 6
  • 243
  • 355
Craig Carter
  • 4,416
  • 16
  • 28
  • Looks similar to the question I asked at https://mathematica.stackexchange.com/questions/112887/ndsolve-mixing-many-scalar-and-vector-equations, which I never did pursue further, so I'm interested in the answer. – SPPearce Oct 23 '18 at 19:56

1 Answers1

4

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

enter image description here

Carl Woll
  • 130,679
  • 6
  • 243
  • 355
  • 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