Error with ode45 when solution approaches 0

Hi guys, I am trying to create an orbital model, but I get this warning when my orbiting particle approaches the zero position: 'Warning: Failure at t=2.334452e+02. Unable to meet integration tolerances without reducing the step size below the smallest value allowed (4.547474e-13) at time t'. I get the same error on the y coordinate and I'm not able to plot a nice elliptic orbit.

G = 6.7e-11;
Ito_mass=3.5e10;
Sat_mass=1.33;

P0=[60,20,0];
a=norm(P0);
v0=sqrt(G*Sat_mass*Ito_mass/a);
vy=v0/sqrt(1+P0(2)/P0(1));
vx=-vy*P0(2)/P0(1);
V0=[vx,vy,0];
t_per=[1:1:700];

initials=[P0(1),V0(1)];
options = odeset('RelTol',1e-4,'AbsTol',1e-6);
[tx,x]=ode45(@myode,t_per,initials,options);
figure
plot(tx,x)

initials=[P0(2),V0(2)];
[ty,y]=ode45(@myodey,t_per,initials,options);
figure
plot(ty,y)

and this is the function I pass to ode45

function myode=myode(t,y)
    G = 6.7e-11;
    Ito_mass=3.5e10;    
    Sat_mass=1.33;
    myode(1)=y(2);
    myode(2)=-G*Ito_mass*Sat_mass*y(1)/(norm(y(1)))^3;
    myode=[myode(1),myode(2)]';

Answer 1

[James Tursa]

You've got a 3D 2nd order problem defined (orbit in 3D space), but you are calling ode45 as if you only had a 1D problem to solve. Also you have a function named the same as the return variable, myode. So pass the entire 3-element position and 3-element velocity (6 elements total) in as initial conditions, and then in your derivative routine calculate the entire 6-element derivative and fix up the names. E.g.,

G = 6.7e-11;
Ito_mass = 3.5e10;
Sat_mass = 1.33;

P0 = [60,20,0];
a = norm(P0);
v0 = sqrt(G*Sat_mass*Ito_mass/a);
vy = v0/sqrt(1+P0(2)/P0(1));
vx = -vy*P0(2)/P0(1);
V0 = [vx,vy,0];
t_per = 1:1420; % <-- Increased the total time

initials=[P0,V0]; % <-- All 6 initial conditions
options = odeset('RelTol',1e-4,'AbsTol',1e-6);
[tx,x] = ode45(@myode,t_per,initials,options);
figure
plot(x(:,1),x(:,2)) % First column is X position, 2nd column is Y position
daspect([1 1 1]); % Set aspect to be the same on all axes

And the deriv routine:

function dy=myode(t,y) % <-- Change the return variable name to dy
    G = 6.7e-11;
    Ito_mass = 3.5e10;    
    Sat_mass = 1.33;
    dy = zeros(size(y));
    dy(1:3) = y(4:6); % Position derivative is simply the velocity
    dy(4:6) = -G * Ito_mass * Sat_mass * y(1:3)/(norm(y(1:3)))^3; % Velocity derivative is gravity

[Dario Lusso]

Thank you James! I have a question tho, I tried to increase the total time to 14200 to have more than one orbit and the following orbits are slightly larger. Do you know why is that? I think it may be due to the solver approximation, but I want to be sure it's not because of a mistake in my formulas. Thanks again!

[James Tursa]

Solver integration error will build up over time, so results that differ from pure constant ellipse will occur.

[Dario Lusso]

Ok, Thanks! One last question, I tried playing around with different P0 values and I get funny results for [20,40,0].. I really can't understand why the particle slows down to almost zero velocity and doesn't run the full orbit.

 %DATA & Initial Conditions
  G = 6.7e-11;
  Ito_mass = 3.5e10;
  Sat_mass = 1.33;
  P0 = [20,40,0];
  a = norm(P0);
  v0 = sqrt(G*Sat_mass*Ito_mass/a);
  vy = v0/sqrt(1+P0(2)/P0(1));
  vx = -vy*P0(2)/P0(1);

if P0(1)==0
    vy=0;
    vx=v0;
end
if P0(2)==0
    vx=0;
    vy=v0;
end
V0 = [vx,vy,0];

t_per = 1:3000; % <-- Increased the total time
initials=[P0,V0]; % <-- All 6 initial conditions

%System Solution
options = odeset('RelTol',1e-4,'AbsTol',1e-6);
[tx,x] = ode45(@myode_corr,t_per,initials,options);
figure
plot(x(:,1),x(:,2)) % First column is X position, 2nd column is Y position
daspect([1 1 1]); % Set aspect to be the same on all axes
figure
plot3(x(:,1),x(:,2),x(:,3));
hold on
[xs,ys,zs]=sphere(100);
surf(xs,ys,zs);
axis([-100 100 -100 100 -100 100])
figure
subplot(4,1,1);
plot(tx, x(:,1));
xlabel('Time [s]', 'FontSize', 10);
ylabel('x [m]', 'FontSize', 10);
subplot(4,1,2);
plot(tx,x(:,2));
xlabel('Time[s]', 'FontSize', 10);
ylabel('y [m]', 'FontSize', 10);
subplot(4,1,3);
plot(tx, x(:,4));
xlabel('Time[s]', 'FontSize', 10);
ylabel('Vx [m/s]', 'FontSize', 10);
subplot(4,1,4);
plot(tx,x(:,5));
xlabel('Time[s]', 'FontSize', 10);
ylabel('Vy [m/s]', 'FontSize', 10);

[James Tursa]

This latest case has a larger initial velocity than the first case. Even though the initial position is smaller, it appears to have an orbit with a longer period because of this larger initial velocity. Orbits in general will conserve energy, so as the particle get farther away from the attracting body it increases potential energy and decreases kinetic energy. This is to be expected. If you run it longer, as long as you haven't given the particle an escape velocity it will eventually start to fall back towards the attracting body. When that happens potential energy will decrease and kinetic energy will increase.