Integration of system of ODEs fails, unable to meet integration tolerances
2 vues (au cours des 30 derniers jours)
Afficher commentaires plus anciens
type556R
le 1 Sep 2022
Commenté : Bjorn Gustavsson
le 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.
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/1114095/image.png)
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 commentaire
Bjorn Gustavsson
le 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.
Réponse acceptée
Bjorn Gustavsson
le 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 commentaires
Bjorn Gustavsson
le 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.
Plus de réponses (1)
Sam Chak
le 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)
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/1114505/image.png)
governs the motion of point mass
(Rocket, in your case) with respect to point mass
(Earth), only for a short period of time.
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/1114510/image.png)
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/1114515/image.png)
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
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/1114520/image.png)
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/1114525/image.png)
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/1114530/image.png)
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.
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/1114535/image.png)
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/1114540/image.png)
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/1114545/image.png)
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
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/1114550/image.png)
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/1114555/image.png)
3 commentaires
Sam Chak
le 2 Sep 2022
Thanks for your clarification. Have you found the reason why it still falls to Earth despite it escapes to the orbit?
Voir également
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!