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.her...@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.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
