Integration of system of ODEs fails, unable to meet integration tolerances
3 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.
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.
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