ode45 orbit glitch

Good morning

I am trying to create a script to integrate the equation r'' = -mu/r^3, by using ode45 and get an elliptical orbit but something is wrong. The orbit is not closing. I believe the initial conditions are not being sent to the function therefore it can't converge.

I call it this way

if true
  % %%
tstep = [0; 100*43000]; %periods
x10 = 10016.34;
x20 = -17012.52;
x30 = 7899.28;
x40 = 2.50;
x50 = -1.05;
x60 = 3.88;
x0 = [x10 x20 x30 x40 x50 x60]';
[t,x] = ode45(@integrator,tstep,x0*1e+3);
plot3 (x(:,1),x(:,2),x(:,3));
grid
end

And the function is

if true
  % function dr_dt = integrator(t,r)
mu = 3.986004e+14;
x = r(1);
y = r(2);
z = r(3);
u = r(4);
v = r(5);
w = r(6);
dx1dt = u;
dx2dt = v;
dx3dt = w;
dv1dt = -mu*x/(((x^2+y^2+z^2)^0.5)^3);
dv2dt = -mu*y/(((x^2+y^2+z^2)^0.5)^3);
dv3dt = -mu*z/(((x^2+y^2+z^2)^0.5)^3);
dr_dt = [dx1dt;%x'
  dx2dt;%y'
  dx3dt;%z'
  dv1dt;%u'
  dv2dt;%v'
  dv3dt];%w'

end

So terra is the orbit file and integrator is the function file. I don't know what i'm doing wrong. I'd appreciate some help.

Alex

Answer 1

[Jan]

 I believe the initial conditions are not being sent

Fortunately you do not have to believe anything in Matlab, but you can use the debugger to check this. But even the output of the integrator shows directly, that the trajectories start at the wanted initial values.

The descriptions "orbit is not closing" and "it can't converge" are not clear. Therefore I do not know, what you assume to be "wrong". When I run your code I get a warning:

 Warning: Failure at t=4.027359e+006.  Unable to meet integration tolerances
 without reducing the step size below the smallest value allowed (7.450581e-009)
 at time t.

The trajectory looks much nicer if you use a smaller tolerance:

opts  = odeset('reltol', 1e-6, 'abstol', 1e-6);
[t,x] = ode45(@integrator,tstep,x0*1e+3, opts);

You cannot expect a numerical integrator to calculate the exact results. You need to adjust the tolerances to you needs.

[Walter Roberson]

The failure to meet integration tolerances suggests you are encountering a singularity. Your equations are singular if the trajectory passes through (0, 0, 0)

[Alexandre Alves Silva]

It's worked, at least until the 100 orbits. It should "close" (instead of falling in the "Earth") because there's no perturbation, but it really was a numerical problem. I just learned this odeset stuff, thanks!