Integration of system of ODEs fails, unable to meet integration tolerances
2 views (last 30 days)
Show older comments
type556R
on 1 Sep 2022
Commented: Bjorn Gustavsson
on 2 Sep 2022
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;
end
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
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.
Accepted Answer
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.
HTH
3 Comments
Bjorn Gustavsson
on 2 Sep 2022
Good. But that still seems to be into the ground. Check:
e_v0 = [x0(3:4);0]/norm(x0(3:4));
e_r0 = [x0(1:2);0]/norm(x0(1:2));
theta_launch = acos(dot(e_r0,e_v0))
That has to be smaller than pi/2.
More Answers (1)
Sam Chak
on 2 Sep 2022
Hi @Luca Biddau
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


3 Comments
Sam Chak
on 2 Sep 2022
Thanks for your clarification. Have you found the reason why it still falls to Earth despite it escapes to the orbit?
See Also
Categories
Find more on Ordinary Differential Equations in Help Center and File Exchange
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!