Can RecurrenceTable make use of CompiledFunction?

Question

I'm trying to implement a discrete-time 2D Verlet algorithm for a point-mass subject to a softened gravitational interaction as a test for a more computationally intensive simulation using RecurrenceTable to recursively compute

$$ \vec{x}_{n+1} = 2 \vec{x}_{n} - \vec{x}_{n-1} + \vec{a}_{n}dt^2 $$

Initial conditions are pseudorandomly generated in a rather complicated way that is not really relevant to the issue I'm experiencing, so consider for instance

Sinitialposition = Developer`ToPackedArray[{-5.35386, -73.709}];
preSinitialposition = Developer`ToPackedArray[{-5.45406, -73.8311}];

Sample values for constants used in the following are

dt = 1.*^-3;
G = 4.49*^3;
M = 1.;
S = 1.;
ϵ = 2.;

Acceleration is given by the CompiledFunction object SAcceleration, defined by

With[{G = G, M = M, S = S, ϵ2 = ϵ^2},
 SAcceleration = 
   Compile[{{SPosition, _Real, 
      1}}, (-G (M + S))/(SPosition.SPosition + ϵ2)^(3/2) SPosition]];

(I have been suggested this syntax in order to avoid MainEvaluate calls while still having the freedom to call for global variables)

When I call

RecurrenceTable[{SR[n + 1] == 
   2 SR[n] - SR[n - 1] + dt^2 SAcceleration[SR[n]], 
  SR[0] == Sinitialposition, SR[-1] == preSinitialposition}, SR, {n, 
  1, 1000}]

I keep running into a warning which tells me, I think, that Mathematica is using the non-compiled version:

CompiledFunction::cfta: "Argument SR[n] at position 1 should be a rank 1 tensor of \!\(\"machine-size real number\"\)s."

I tried investigating using Trace, and the relevant portion of the output (at position [[1, 1, 2, 2, 2, 2]]) is:

CompiledFunction[{SPosition},-(((4490. 2.) SPosition)/(SPosition.SPosition+4.)^(3/2)),-CompiledCode-][SR[n]]

So it seems that RecurrenceTable calls the function with the given symbolic arguments. Does this mean I cannot make use of CompiledFunction objects or does RecurrenceTable auto-compile? Am I going to encounter performance hits because of this? Are there any workarounds?

Answer 1

The reason for this message is that the compiled function is called with the symbolic argument SR[n] in the definition of the recurrence relation:

SAcceleration[SR[n]]

CompiledFunction::cfta: "Argument SR[n] at position 1 should be a rank 1 tensor of machine-size real numbers."

-((8980. SR[n])/(4. + SR[n].SR[n])^(3/2))

The recurrence is then evaluated using the symbolic form returned by SAcceleration. While I don't think you'll gain much (if anything) in the way of performance by using a compiled function for such a simple expression, if you want to get rid of the message and avoid conversion to the uncompiled form, you can prevent SAcceleration from evaluating for symbolic arguments. One way to do this is shown by @Sjoerd, but in version 8 there is another way that's specific to compiled functions, using RuntimeOptions as follows:

RuntimeOptions -> {"EvaluateSymbolically" -> False}

Adding this we have:

With[{G = G, M = M, S = S, ϵ2 = ϵ^2},
 SAcceleration = Compile[{{SPosition, _Real, 1}},
  (-G (M + S))/(SPosition.SPosition + ϵ2)^(3/2) SPosition,
  RuntimeOptions -> {"EvaluateSymbolically" -> False}
 ]
];

Now, SAcceleration[SR[n]] returns unevaluated and the resulting behaviour is as you intended.

Answer 2

If you could define your compiled function using an intermediate step I think it should work:

With[{G = G, M = M, S = S, ϵ2 = ϵ^2}, 
  SAcceleration1 = 
      Compile[{{SPosition, _Real,1}}, 
              (-G (M + S))/(SPosition.SPosition + ϵ2)^(3/2) SPosition]
];

SAcceleration[{x_?NumericQ, y_?NumericQ}] := SAcceleration1[{x, y}]

Running it:

rt = RecurrenceTable[{SR[n + 1] == 2 SR[n] - SR[n - 1] + dt^2 SAcceleration[SR[n]], 
                      SR[0] == Sinitialposition, 
                      SR[-1] == preSinitialposition}, 
                      SR, {n, 1, 1000}
     ];

ListPlot[Norm /@ Differences@Partition[rt, 2, 1]]