two body problem using ode45

Question

siyeong Jang 2019-11-19

1
链接

此问题的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/491933-two-body-problem-using-ode45

回答： Meysam Mahooti 2021-7-28

% Assumption that earth is orbiting circulary around sun

v=0.0172; % angular velocity of earth [au/day]

p=365; % orbital time of earth [day]

a=1; % orbital radius of earth [au]

[t,x]=ode45(@earth,[0:0.001:p],[a;0;0;v]);

plot(x(:,1),x(:,3));

function dxdt=earth(t,x)

g=2.959*10^(-4);

m=1;

dxdt=[x(2);-g*m*x(1)/(x(1)^2+x(3)^2)^1.5;x(4);-g*m*x(3)/(x(1)^2+x(3)^2)^1.5];

end

theoretically, the earth must moves closed orbit.

but it doesn`t for my code

long time later the earth is converge into sun

what does matter?

0 个评论
显示 -2更早的评论隐藏 -2更早的评论

请先登录，再进行评论。

请先登录，再回答此问题。

Answer 1

James Tursa 2019-11-19

1
链接

此回答的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/491933-two-body-problem-using-ode45#answer_402275

编辑：James Tursa 2019-11-19

在 MATLAB Online 中打开

The RK methods don't match up well with the orbit DE problem because the integration errors tend to be systematic (e.g., always integrates a bit low) and build up over time. You can counteract this by specifying a tighter tolerance. Also, to get the best precision for your initial circular orbit calculate v programatically instead of hard coding it. E.g.,

g = 2.959e-4;
v = sqrt(g/a);
n = sqrt(g/a^3);
p = 2*pi/n;
options = odeset('RelTol',1e-10);
[t,x]=ode45(@earth,[0:0.001:20*p],[a;0;0;v],options);

2 个评论
显示无隐藏无

siyeong Jang 2019-11-20

Thank you for your advice.

I have one more question.

what is the diffrent between solver ode45 and ode113.

it`s make different result in restricted three body problem.

whait is the correct? 제목 없음.png

James Tursa 2019-11-20

I could only point you to the doc for details of the differences:

https://www.mathworks.com/help/matlab/ref/ode113.html?searchHighlight=ode113&s_tid=doc_srchtitle

In particular, see the "ODE with Stringent Error Threshold" example for a comparison of ode45 and ode113 in a two-body problem.

请先登录，再进行评论。

Answer 2

Jim Riggs 2019-11-19

1
链接

此回答的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/491933-two-body-problem-using-ode45#answer_402295

编辑：Jim Riggs 2019-11-19

I have studied the 2-body problem and the accuracy of different numerical solutions. See my comments in this thread:

Discussion on Euler's Method - 2 body problem example

I have found that the 4th order runge kutta is the most efficient solver for the 2 body problem. (This means that ode45 is a good choice)

I will reitterate what James Tursa says, you will NEVER get a closed orbit using a numerical approximation, it is only a matter of controlling the size of the error. This is why we do sensitivity studies to measure the numerical error as a functionof solver type and time step. The 2-body problem is a good one for evaluating numerical erors, since the 2 body problem SHOULD be a closed orbit, any deviation is ascribed to the solver, (but, note this is only true when positions are mesured from the system barycenter. Since you have pinned the sun position so that it does not move, you are OK)

Again, see the discussion at the given link.