Numerical representation

[Jon Herman]

And for the sake of additional completeness (I'm sorry I didn't think of all
this in one go): my derivative function in Python produces results that
agree with MATLAB to order e-16 (machine precision), so the error is
definitely building up in my integrator.



On Mon, Mar 7, 2011 at 11:59 AM, Jon Herman <jfc.herman at gmail.com> wrote:

> 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/79103399/attachment.html>