Take the 2-minute tour ×
Mathematica Stack Exchange is a question and answer site for users of Mathematica. It's 100% free, no registration required.

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.;
\[Epsilon] = 2.;

Acceleration is given by the CompiledFunction object SAcceleration, defined by

With[{G = G, M = M, S = S, \[Epsilon]2 = \[Epsilon]^2},
 SAcceleration = 
   Compile[{{SPosition, _Real, 
      1}}, (-G (M + S))/(SPosition.SPosition + \[Epsilon]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?

share|improve this question
    
If RecurrenceTable[] were only HoldAll... too bad it isn't, so there you are. –  J. M. Jul 1 '12 at 12:54
    
@OleksandrR. It works. Why don't you state that as an answer? It is precisely what I've been looking for! :O –  Andrea Colonna Jul 1 '12 at 14:40
add comment

2 Answers

up vote 11 down vote accepted

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, \[Epsilon]2 = \[Epsilon]^2},
 SAcceleration = Compile[{{SPosition, _Real, 1}},
  (-G (M + S))/(SPosition.SPosition + \[Epsilon]2)^(3/2) SPosition,
  RuntimeOptions -> {"EvaluateSymbolically" -> False}
 ]
];

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

share|improve this answer
add comment

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

With[{G = G, M = M, S = S, \[Epsilon]2 = \[Epsilon]^2}, 
  SAcceleration1 = 
      Compile[{{SPosition, _Real,1}}, 
              (-G (M + S))/(SPosition.SPosition + \[Epsilon]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]]

Mathematica graphics

share|improve this answer
    
+1. I'm too late, having some chat when editing, LOL –  Silvia Jul 1 '12 at 13:27
    
@Silvia Happens all the time. If I hadn't been editing the question first I'd be quicker too. I'll vote for you as well, as your answer is equally valid. –  Sjoerd C. de Vries Jul 1 '12 at 13:29
    
It is indeed equipollence. Maybe I delete mine. (if it won't do harm for others~) –  Silvia Jul 1 '12 at 13:34
    
@Silvia The main difference would be that your solution would fail in the unlikely case that someone would enter a list consisting of symbolic values. –  Sjoerd C. de Vries Jul 1 '12 at 13:38
    
Yes, I got lazy for it :D And NumericQ would not be enough for preventing a Complex list :) Maybe something like x_?(Element[#, Reals] || Element[#, Integers] &) –  Silvia Jul 1 '12 at 13:48
show 3 more comments

Your Answer

 
discard

By posting your answer, you agree to the privacy policy and terms of service.

Not the answer you're looking for? Browse other questions tagged or ask your own question.