Integration of system of ODEs fails, unable to meet integration tolerances

2 views (last 30 days)
I'm trying to integrate the following system of four differential equations related to the two body problem, for an escape trajectory from Earth:
tspan = [0 9.944e+04]; %[s]
x0 = [-5826; -3249; 9.887; 5.51]; %position[km] and speed[km/s] in ECI frame after burn out
[t,x] = ode15s(@ODE_departure,tspan,x0);
function dxdt = ODE_departure(tspan,x) %two body problem diff. equations
mu = 398600; %Earth gravitational parameter in km3/s2
dxdt = zeros(4,1);
r = x(1:2);
dxdt(1:2) = x(3:4);
dxdt(3:4) = -mu*x(1:2)/norm(r)^3;
What I end up is:
  • The error "Warning: Failure at t=3.957528e+02. Unable to meet integration tolerances without reducing the step size below the smallest value allowed (9.094947e-13) at time t"
  • A trajectory that "goes to zero", to the center of the Earth, meaning the solutions x(:,1) and x(:,2) gets smaller and smaller when they should grow bigger instead
  • Speeds, so solution x(:,3),x(:,4) that constantly grow to 10^5 km/s, when they should instead converge to about 2km/s
I already tried using the more classic ode45, or reducing the tolerances, but basically nothing changes.
I thought the problem was how I re-wrote the second order differential equation , but I checked it and wrote it in different ways many times.
I think that the integration is failing because when the distance r (norm(x(:,[1,2])))goes to zero, there will be a singularity. But I don't understand why this is happening during the integration.
This is my first question, I hope it's clear enough, if needed I'll provide graphical results and more informations, thank you.
  1 Comment
Bjorn Gustavsson
Bjorn Gustavsson on 2 Sep 2022
To check that you get away in the right direction, limit the integration to a short enough period of time that the solver doesn't get itself into trouble. Then plot the initial trajectory and the sphere of Earth, "rocket" should head out, right? And definitely not intersect the surface. You should make checking these things an instinctive/habitual behaviour, also check the conservation of energy of the rocket in this case, or other known constants of motion (or the similar), when using numerical integration-schemes.

Sign in to comment.

Accepted Answer

Bjorn Gustavsson
Bjorn Gustavsson on 1 Sep 2022
After only quickly looking at your problem, it seems to me that you launch the body in a direction into Earth. The first 2 and last 2 components of x0 are pretty antiparallel. Try to fix that - either by moving the launch-point or change the direction of the launch-velocity:
x0 = [-5826; -3249; -9.887; -5.51]; % for example...
I think this should be enough to solve your issue.
type556R on 2 Sep 2022
I solved it, thank you. I corrected that again by first computing the right velocity direction. Another cause of the problem was writing
dxdt(3:4) = -mu*x(1:2)/sqrt(x(1)^2+x(2)^2)^3;
Instead of the correct line:
dxdt(3:4) = -mu*x(1:2)/norm([x(1:2)])^3;
I have no idea why that was making the trajectory go away from Earth and then straight to its center, I think should be equivalent to norm([x1,x2])

Sign in to comment.

More Answers (1)

Sam Chak
Sam Chak on 2 Sep 2022
There is nothing wrong with the assumption of the equation. However, the idealized differential equation of the Keplerian two-body motion (in acceleration mode)
governs the motion of point mass (Rocket, in your case) with respect to point mass (Earth), only for a short period of time.
If you run the simulation, it will always "fall back" to the center of Earth because the equation lacks the perturbative acceleration acting upon (Rocket). So, the true equation should be
What is ? It can be any perturbing effect that induces the acceleration. For example, the non-sphericity of the primary mass (aka, the effect), the presence of the gravitational fields from other bodies (Sun and Moon), atmospheric drag, etc.
But in your case, the most significant component of comes from your Rocket, the Rocket engine should continuously burn and expel propellant to provide acceleration to the rocket to counter Earth's gravity until achieving the Escape Velocity (and maintaining its speed until it reaches the orbit). If other insignificant perturbing effects are ignored, then the equation becomes
type556R on 2 Sep 2022
Yes, the initial velocity vector was oriented in the wrong direction, and there was a wrong line in the definition of the ODEs to be integrated, you can find the details in the comments of the accepted answer if you're interested.

Sign in to comment.

Community Treasure Hunt

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

Start Hunting!