Coupled second order ODE with three variables

1 vue (au cours des 30 derniers jours)
Erfan Saeedi
Erfan Saeedi le 27 Jan 2020
Commenté : Erfan Saeedi le 28 Jan 2020
Hi, I am trying to solve this system of equations (above).
I have writen written these codes for it but i'm not getting correct answers.
r0=[1000 0 1000 0 1000 0];
[t,r]=ode45(@f,[0 tmax],r0);
function dr=f(~,r)
G=6.6742*(10^(-11));
M_e=5.972*(10^24);
%r(1)=x
%r(3)=y
%r(5)=z
dr = zeros(6,1);
dr(1)=r(2);
dr(2)=-(G*M_e)*r(1)/(r(1)^2+r(3)^2+r(5)^2)^(3/2);
dr(3)=r(4);
dr(4)=-(G*M_e)*r(1)/(r(1)^2+r(3)^2+r(5)^2)^(3/2);
dr(5)=r(6);
dr(6)=-(G*M_e)*r(5)/(r(1)^2+r(3)^2+r(5)^2)^(3/2);
end
I'd very much appreiciate it if you could help me with this problem.

Réponses (1)

James Tursa
James Tursa le 28 Jan 2020
Modifié(e) : James Tursa le 28 Jan 2020
Problems:
1) Units. You really should annotate ALL of your constants with descriptions and units so that the reader knows what these magic numbers are. E.g.,
G = 6.6742*(10^(-11)); % Gravitational Constant, ( m^3 / ( kg * s^2) )
M_e = 5.972*(10^24); % Mass of the Earth, (kg)
2) Using those constants for G and M_e with those units, this dictates what the units of r must be in order for your derivative equations to work. r must be (m) for position and (m/s) for velocity. That means your initial starting position of [100, 100, 100] is about 173 meters from the center of the Earth ... deep inside the nickel/iron core! Plus, this doesn't even match your picture with starting position of [1000, 1000, 1000]. So I'm not sure what is really intended here.
3) And it has a 0 starting velocity. This makes no sense for an orbit problem ... it would be a rectilinear orbit and fall straight down into the Earth from the start.
I'm guessing that the 1000 for position is probably meant to be an altitude. Then you need to add a realistic velocity. For a circular orbit we can just go with v = sqrt(mu/r).
Summary of the corrections would be something like this up front in your code::
G = 6.6742*(10^(-11)); % Gravitational Constant, ( m^3 / ( kg * s^2) )
M_e = 5.972*(10^24); % Mass of the Earth, (kg)
R_e = 6378137; % Equatorial radius of Earth, (m)
mu = G * M_e; % mu of the Earth (m^3 / s^2)
h = 1000; % Initial altitude of each axis, (km)
rx = (h*1000 + R_e); % Radius of each axis, (m)
r = sqrt(3) * rx; % Total radius for an [rx,rx,rx] position, (m)
vcirc = sqrt(mu/r); % Circular velocity at initial point, (m/s)
vx = vcirc/sqrt(2); % Circular velocity distributed for two axes (m/s)
r0 = [rx vx rx -vx rx 0]; % Initial states for circular orbit, (m) and (m/s)
[t,r]=ode45(@f,[0 tmax],r0);
I just picked a somewhat arbitrary starting position since I'm not sure what you really want for that [1000, 1000, 1000] that you list. I also picked an arbitrary direction for the velocity that was perpendicular to the radius vector in order to get a circular orbit.
SIDE NOTE: In the future, it is usually easier to group all the position elements next to each other and all of the velocity elements next to each other. It is easier to vectorize the equations & derivatives that way.
  3 commentaires
James Tursa
James Tursa le 28 Jan 2020
I seriously doubt that the initial velocity is really 0. You have in all likelihood misunderstood the initial conditions of the problem.
Erfan Saeedi
Erfan Saeedi le 28 Jan 2020
the problem isn't realistic. and i didn't choose the initial conditions.

Connectez-vous pour commenter.

Catégories

En savoir plus sur Programming dans Help Center et File Exchange

Community Treasure Hunt

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

Start Hunting!

Translated by