Re: MATLAB 5.3.1 (1999) and latest Octave math discrepancy ..

[Stephen Montgomery-Smith]

First, it is very likely that the equation you are solving is sensitive
to floating point errors.  That is to say, a small numerical difference
(maybe in the last bit of the mantissa) might later on be greatly
exaggerated.  This is a well known phenomenon in solving ordinary
differential equations (ODE), and it would be hard to say which, if any,
is giving the correct answer - MATLAB or octave.  Another equation which
illustrates this more profoundly is the famous Lorenz attractor
http://en.wikipedia.org/wiki/Lorenz_system

Secondly, your method of solving the (ODE) appears to be the so called
Euler Method http://en.wikipedia.org/wiki/Euler_method.  It is well
known to be rather inaccurate (hence your need for a very small time
step).  Both octave and matlab have much more sophisticated ODE solvers,
called lsode and ode45 respectively.  Or if you want to write your own
solver, I would suggest the Runge-Kutta Method of order 4 with a fixed
step size http://mathworld.wolfram.com/Runge-KuttaMethod.html.  It is
actually fairly straightforward to program, and it will give far more
accurate answers than your program.



On 03/22/2013 07:27 PM, Will Flannery wrote:
> I have no problem showing the code.   First I'll describe what it does.
> 
> The program simulates the Juno space probe ... the probe is launched from a
> low earth satellite orbit and goes into a 2 year orbit around the sun, after
> 2 years it returns to earth, trailing it by just a bit, and earth's
> gravitational field pulls it, increasing it's velocity so that it breaks its
> 2 year orbit and is 'slingshot' toward Jupiter.  You can see a simulation
> here ...
> 
> http://www.youtube.com/watch?v=sYp5p2oL51g
> 
> The program computes the gravitational forces of the earth and the sun on
> the probe, resolves them into x and y components, an projects them over 2
> second subintervals, it does the same thing for the earth's orbit around the
> sun.
> 
> I did the program a long time ago using MATLAB, and just downloaded
> Octave3.6.2gcc4.6.2 a few days ago.  I cut and pasted the program into
> Octave and it ran no problem.
> 
> The thing looks good until the probe returns to earth, and there the earth's
> gravity takes over, and the differences in the graphs, unnoticeable to this
> point, get magnified greatly by earth's gravity and the trajectories diverge
> significantly.
> 
> The code - for some reason way back when, when I was trying to get this
> thing to work, I broke the code into two scripts, which I have preserved
> just because I don't want to change anything.  I ran it with MATLAB and got
> the result shown.
> 
> Note - the program can't store every state value for graphing, as it
> requires too much storate, so it stores the state every 2000 iterations. 
> This produces a lot of unessential code.
> 
> The main program ... 
> _______________________  Orbit11  ____________________________
> 
> clear
> 
> 
> G=6.7e-11;
> mSun = 1.9e30;   
> rEarthOrbit = 150e9;
> rSun=695500000;
> mEarth = 5.9742e24;   
> rEarth = 6.378e6;
> 
> dt=2;
> N=3*365*60*60*24 / dt
> 
> M=2000;
> I=N/M+1;
> Ii=1;
> 
> XE=zeros(I,1); YE=zeros(I,1); %VXE=zeros(I,1); VYE=zeros(I,1);
> XP=zeros(I,1); YP=zeros(I,1); %VYP=zeros(I,1); VYP=zeros(I,1); 
> RPE=zeros(I,1);VP=zeros(I,1);
> rPEmin=200000;
> rPE(1)=rPEmin;
> 
> t1=0; T(1)=0;
> 
> % set launch parameters for a earth orbit
> 
> xE1=rEarthOrbit; XE(1)=xE1;         % starting x coordinate
> vxE1=0; vXE(1)=vxE1;            % starting x velocity
> yE1=0;  yE(1)=yE1;          %2*rEarthOrbit*pi/(365*60*60*24)
> vyE1=28500;  vYE(1)=vyE1;               % starting y velocity
> 
> % set launch parameters for probe 
> xP1=rEarthOrbit + rEarth + 200000;  XP(1)=xP1;        
> vxP1=0;   vXP(1)=vxP1;       
> yP1=0;  yP(1)=yP1;
> LV=5265;  % 5260 gives slingshot with collision  
> 
> vyP1=28500 + LV;  vYP(1)=vyP1; vP(1)=vyP1;   % earth's velocity + launch
> velocity
> 
> Orbit11b
> 
> length(XP)
> I2=30400;
> I1=30100;
> 
> 
> plot(XP,YP,'r')                         % plot the probe's orbit
> 
> hold on
> plot(XE,YE)                         % plot the earth's orbit
> 
> __________ Orbit11b __________________
> 
> for i=1:N
>    % project the earth's orbit
> 
>    t2 = t1+ dt;
>    xE2=xE1+vxE1*dt;        % project x position
>    yE2=yE1+vyE1*dt;        % project y position
>    
>    % resolve the earth sun gravity vector
>    rE = sqrt(xE1^2 + yE1^2); % compute rE
>    g = -G*mSun/rE^2;             % compute gravity magnitude
>    ctheta = xE1/rE;             % compute cos theta
>    stheta = yE1/rE;             % compute sin theta
>    
>    % now we're ready to project velocity
>    vxE2= vxE1+ g*ctheta*dt;  % project x velocity
>    vyE2 = vyE1+ g*stheta*dt;  % project y velocity
>    
>    % project the probe's orbit
> 
>    xP2=xP1+vxP1*dt;        % project x position
>    yP2=yP1+vyP1*dt;        % project y position
>    
>    % resolve the probe sun gravity vector
>    rP = sqrt(xP1^2 + yP1^2); % compute rP
>    g = -G*mSun/rP^2;             % compute gravity magnitude
>    ctheta = xP1/rP;             % compute cos theta
>    stheta = yP1/rP;             % compute sin theta
>    
>    % resolve the probe earth gravity vector
>    xPE = xP1 - xE1;
>    yPE = yP1 - yE1;
>    rPE = sqrt((xPE^2 + yPE)^2); % compute rPE
>    gE = -G*mEarth/rPE^2;             % compute gravity magnitude
>    comega = xPE/rPE;             % compute cos theta
>    somega = yPE/rPE;             % compute sin theta
> 
>    
>    % now we're ready to project velocity
>    vxP2 = vxP1+ (g*ctheta + gE*comega)*dt;  % project x velocity
>    vyP2 = vyP1+ (g*stheta + gE*somega)*dt;  % project y velocity
>       
>    if mod(i,M)== 0        
>       Ii=Ii+1;
>  
>       Ii*M*dt/(60*60*24)
>       XP(Ii)= xP2;
>       %VXP(Ii)= vxP2;
>       YP(Ii)= yP2;
>       %VYP(Ii)= vyP2;
>       
>       XE(Ii)= xE2;
>       %VXE(Ii)= vxE2;
>       YE(Ii)= yE2;
>       %VYE(Ii)= vyE2;
>       
>       VP(Ii)=sqrt(vxP2^2+vyP2^2);
>       RPE(Ii)=rPEmin;
>       rPEmin=rPE;
>       end
>    
>    if(rPE<rPEmin)rPEmin=rPE;
>    end
>    
>   
> 
>    t1=t2;
>    xP1=xP2;
>    vxP1=vxP2;
>    yP1=yP2;
>    vyP1=vyP2;
>    xE1=xE2;
>    vxE1=vxE2;
>    yE1=yE2;
>    vyE1=vyE2;
>   _____________________________________ 
> 
> 
>    
> 
> 
> 
> 
> 
> 
> 
> --
> View this message in context: 
> http://octave.1599824.n4.nabble.com/MATLAB-5-3-1-1999-and-latest-Octave-math-discrepancy-tp4651123p4651130.html
> Sent from the Octave - Maintainers mailing list archive at Nabble.com.
> 
>