Numerical representation

Jon Herman jfc.herman at gmail.com
Mon Mar 7 13:59:35 EST 2011


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.herman at 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!
>
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://mail.python.org/pipermail/python-list/attachments/20110307/68acbdfb/attachment-0001.html>


More information about the Python-list mailing list