And for the sake of completeness, the derivative function I am calling from
my integrator (this is the 3 body problem in astrodynamics):
def F(mu, X, ti):
r1= pow((pow(X[0]+mu,2)+pow(X[1],2)+pow(X[2],2)),0.5)
r2= pow((pow(X[0]+mu-1,2)+pow(X[1],2)+pow(X[2],2)),0.5)
Ax= X[0]+2*X[4]-(1-mu)*(X[0]+mu)/r1**3-mu*(X[0]-(1-mu))/r2**3
Ay= X[1]-2*X[3]-(1-mu)*X[1]/r1**3-mu*X[1]/r2**3
Az= -(1-mu)*X[2]/r1**3-mu*X[2]/r2**3
XDelta=array([X[3], X[4], X[5], Ax, Ay, Az])
return XDelta
On Mon, Mar 7, 2011 at 11:50 AM, Jon Herman <jfc.her...@gmail.com> wrote:
> Sorry Robert, I'd missed your post when I just made my last one. The output
> I am getting in Python looks as follows:
>
> array([ 9.91565050e-01, 1.55680112e-05, -1.53258602e-05,
> -5.75847623e-05, -9.64290960e-03, -8.26333458e-08])
>
> This is the final state vector, consisting of 6 states (postion and
> velocity), although this isn't really relevant right now.
>
> In MATLAB, using the same process for the RKF78, I get the following
> output:
>
> Xend =
>
> Columns 1 through 3
>
> 9.915755153307796e-001 3.592556838089922e-016
> -1.523933534321440e-005
>
> Columns 4 through 6
>
> -7.175069559491303e-015 -9.624755221500220e-003
> 1.289789641929434e-017
>
> As you can see, the results are close but there is a big difference in
> precision (the 2nd, 4th and 6th entries are supposed to be zero under the
> intial and final conditions I am running).
> See also the post I just made where you can see the Python code I am using
> (MATLAB code is identical but translated)
>
> This is for a fixed timestep in both Python and Matlab. If I let the code
> run with an identical method for timestep correction (based on a tolerance),
> Python will use a timestep approximately 10^5 SMALLER than Matlab uses. So
> I'm really sure it's not a representation issue, but actually a precision
> issue.
>
> Thanks for any advice!
>
--
http://mail.python.org/mailman/listinfo/python-list
