Two-body problem program gone wrong

46 Ansichten (letzte 30 Tage)
Matthew Lee
Matthew Lee am 24 Feb. 2021
Hi everyone... I have written a simple program on matlab that will draw an orbit using ode45 to integrate the two body problem which is rddot = -mue.*rv/R^3
where:
rddot: second derivative of position vector (acceleration)
mue: Earth constant (which is known to be = 398600 km^3/s^2)
rv: position vector
R: position scalar
and my problem is the orbit will complete one circle without a problem but it will decay as it complete more than one period around the earth which I'm sure it's not the realistic case..
for example, completing one period it will look normal:
if I wanted to complete more than one period, the orbit will decay:
this picture is taken after completing 20 rotations.
someone please help if found any problem in my code or my mathematical expression,
this is my code:
r0=[-7072.7292 3150.61 123.8]'; %km
v0=[0.4600 -4.2446 -7.312]'; %km/s
y0=[r0;v0];
T=14332.25;
tspan=[0 T];
odefun=@dfdf;
[t,y]= ode45(odefun,tspan,y0);
xv = y(:,1);
yv = y(:,2);
zv = y(:,3);
plot3(xv,yv,zv)
axis equal
function [dydt]=dfdf(t,y)
rv=y(1:3);
R=norm(rv);
mue=398600.44189;
rdot=y(4:6);
rddot = -mue*rv/R^3;
dydt = [rdot;rddot];
end
  1 Kommentar
Just Manuel
Just Manuel am 24 Feb. 2021
More a philosophical answer, thus just as a comment:
In my experience, it is very hard to find parameters that produce a stable orbit. The fact that we are orbiting the sun is thus rather a miracle :) In fact, even the distance to the sun is growing (https://www.forbes.com/sites/startswithabang/2019/01/03/earth-is-drifting-away-from-the-sun-and-so-are-all-the-planets).
So I suspect that there is nothing wrong with your simulation. You just have not found the perfect parameter set. Try tweakind them a bit.
Cheers
Manuel

Melden Sie sich an, um zu kommentieren.

Akzeptierte Antwort

Jan
Jan am 24 Feb. 2021
Decrease the tolerance to reduce the truncation errors:
opt = odeset('AbsTol', 1e-8, 'RelTol', 1e-8);
[t,y]= ode45(odefun,tspan,y0, opt);
  1 Kommentar
Matthew Lee
Matthew Lee am 24 Feb. 2021
Thank you very much!!
that worked for me!

Melden Sie sich an, um zu kommentieren.

Weitere Antworten (2)

darova
darova am 24 Feb. 2021
rddot = -mue*rv/R^3;
This line tells that you have non-constant angular acceleration
So it's ok if you orbit is changing

edoardo de angeli
edoardo de angeli am 30 Jul. 2021
How did you manage to iterate the process for more than one orbit?

Kategorien

Mehr zu Programming finden Sie in Help Center und File Exchange

Community Treasure Hunt

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

Start Hunting!

Translated by