[Date Prev][Date Next][Thread Prev][Thread Next][Date Index][Thread Index]
Re: MATLAB 5.3.1 (1999) and latest Octave math discrepancy ..
From: |
Stephen Montgomery-Smith |
Subject: |
Re: MATLAB 5.3.1 (1999) and latest Octave math discrepancy .. |
Date: |
Fri, 22 Mar 2013 23:13:44 -0500 |
User-agent: |
Mozilla/5.0 (X11; Linux i686; rv:17.0) Gecko/20130308 Thunderbird/17.0.4 |
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.
>
>