ode45 orbit glitch

10 次查看（过去 30 天）
Alexandre Alves Silva
编辑： Alexandre Alves Silva 2017-6-6
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

请先登录，再回答此问题。

采纳的回答

Jan
Jan 2017-6-6
编辑：Jan 2017-6-6
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.
  2 个评论
Walter Roberson
Walter Roberson 2017-6-6
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
Alexandre Alves Silva 2017-6-6
编辑：Alexandre Alves Silva 2017-6-6
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!

请先登录，再进行评论。

更多回答（0 个）

类别

在 Help Center 和 File Exchange 中查找有关 Gravitation, Cosmology & Astrophysics 的更多信息

标签

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by