MATLAB: Ode45 orbit glitch

initial conditionsode45orbit

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

Best Answer

 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.
Related Question