d-dimensional compiled function in NDSolve

Question

Based on this answer, I understand how to use a one dimensional compiled function in NDSolve. Mathematica NDSolve and 'Compile'?

I'd like to be able to use a d-dimensional compiled function f, in the same way. f would be a vector-valued function which returns a vector, somthing like this:

f = Compile[{{x, _Real, 1}}, x]

But imagine that it's costly to evaluate.

I've dreamed up some syntax which I know doesn't work, but which kind of says what I want to happen:

rhs[x_?VectorQ[#, NumericQ]] := f[x] 
NDSolve[Join[Thread[Table[x[i]'[t],{i,d}]-rhs[Table[x[i][t],{i,d}]]==0],initialconditions], Table[x[i][t],{i,d}], {t, 0, 1}] 

I know that the above code won't work, because Mathematica won't know how to add together the vector of derivatives and the rhs function until it's been evaluated numerically.

Does anybody know a syntax which works correctly for this?

Note that I can't just break down the vector function f into d separate functions; it's costly to evaluate and I need to make sure that f is called only once each time its components are accessed for a given value of x.

Answer 1

There are several questions on the site that ask about or deal with array-valued variables in differential equations, especially NDSolve/NDSolveValue. They involve other issues as well, and this question puts the basic problem quite simply and clearly. For instance, we can use the "vector-valued solution" from this Q&A as an MWE.

Here's a complete example:

d = 4;  (* dimension *)
rhs[t_, x_?(VectorQ[#, NumericQ] &)] := (Cos[Range[d] t] *
    AdjacencyMatrix@CycleGraph[d, DirectedEdges -> True]).x - x + 1;
sol = NDSolve[{x'[t] == rhs[t, x[t]], x[0] == Range[d]}, x, {t, 0, 2}];

The solution consists of 73 steps, the value at each step a 4-vector:

x["ValuesOnGrid"] /. First@sol // Dimensions
(*  {73, 4}  *)

One constraint of the solution is that evaluating it at, say, t=1.5, results in getting all four components.

x[1.5] /. First[sol]
(*  {1.94518, 1.23931, 1.10914, 1.54423}  *)

So to use the solution takes manipulation. For instance, an easy (if less well-known) way to get a plot is this (from the same Q&A):

ListLinePlot[
 Transpose[{Flatten[x["Grid"] /. sol], #}] & /@
  (Transpose[x["ValuesOnGrid"]] /. First@sol), 
 PlotLegends -> Array[Inactive[Part][x, #] &, x["OutputDimensions"] /. First@sol]]

To get the value of a particular coordinate, one cannot simply execute x[1.5][[1]] /. First@sol, because x[1.5][[1]] is evaluated to 1.5 first: the next step becomes 1.5 /. First@sol, which does not give the value of x. Instead, one can get the order of evaluation right or use Indexed:

(x[1.5] /. First@sol)[[1]]
Indexed[x[1.5], 1] /. First@sol
(*  1.94518  *)

Alternatively, one can reconstruct the individual coordinate functions from sol, which again requires some array gymnastics:

compSol = With[{data = Hold[Transpose[{x["Grid"], ##}] & @@@ 
    Transpose[{x["ValuesOnGrid"], x'["ValuesOnGrid"]}, {2, 3, 1}]] /. First@sol},
  Interpolation /@ ReleaseHold@data]
(* results in a list of 4 InterpolationFunction[]'s *)

ListLinePlot@compSol
(* gives the same graph as above (without the legend) *)

compSol[[1]]@ 1.5  (* evaluate first component at 1.5 *)
(*  1.94518  *)