ODE solver for restricted 3 body problem

Question

Romio 2022-9-19

0
链接

此问题的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/1808140-ode-solver-for-restricted-3-body-problem

回答： Kevin Valencia 2022-12-13

I am trying to implement the following problem(restricted 3 body problem), but I don't know if my implementation is correct since the plot does not look right to me (niether phase plane nor in time). Is my impementation correct? I converted the 2nd order ODE system of 1st order and used the given parameters. Here is the problem and my code

clc,clear,close all
mu = 0.012277471
Y0 = [0.994; 0.0; 0.0; -2.00158510637908252240527862224];
t0 = 0; 
T = 17.0652164601579625588917206249;
tspan = [t0 T];
options = odeset('RelTol',1e-13,'AbsTol',1e-16 * ones(4,1));
[t,Y] = ode15s(@(t,Y) ThreeBodyProblem(t,Y,mu), tspan, Y0,options);
Y = Y';
plot(t,Y)
% plot(Y(2,:),Y(4,:))
function dydt = ThreeBodyProblem(t,y,mu)
dydt = zeros(4,1);
 
mu_p = 1 - mu;
D1 = ( (y(1) + mu)^2 + y(3)^2 )^(3/2);
D2 = ( (y(1) - mu_p)^2 + y(3)^2 )^(3/2);
dydt(1) = y(2);
dydt(2) = y(1) + 2 * dydt(3) - mu_p * ((y(1) + mu)/D1) - mu * ((y(1) - mu)/D2);
dydt(3) = y(4);
dydt(4) = y(3) + 2 * dydt(1) - mu_p * (y(3)/D1) - mu * (y(3)/D2);
end

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

请先登录，再进行评论。

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

Answer 1

James Tursa 2022-9-19

0
链接

此回答的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/1808140-ode-solver-for-restricted-3-body-problem#answer_1056685

编辑：James Tursa 2022-9-19

在 MATLAB Online 中打开

In this line:

dydt(2) = y(1) + 2 * dydt(3) - mu_p * ((y(1) + mu)/D1) - mu * ((y(1) - mu)/D2);

you use dydt(3) before you set its value, so this is incorrect. Simply reorder your equations to avoid this (or use y(4) instead of dydt(3) on rhs). Also it looks like you have some typos with signs and mu prime stuff. E.g.,

dydt(1) = y(2);
dydt(3) = y(4);
dydt(2) = y(1) + 2 * dydt(3) - mu_p * ((y(1) + mu)/D1) - mu * ((y(1) - mu_p)/D2);
dydt(4) = y(3) - 2 * dydt(1) - mu_p * (y(3)/D1) - mu * (y(3)/D2);

Side note: You have 30 digits of precision listed for a couple of your values, but double precision only has about 15 digits of precision, so those trailing digits are meaningless.

2 个评论
显示无隐藏无

Romio 2022-9-19

Oops thanks for catching the typos!! Is there any way to increase the percison to that level (just a side question) in MATLAB?

James Tursa 2022-9-19

编辑：James Tursa 2022-9-19

You would have to set everything up as an extended precision problem, e.g. using the Symbolic Toolbox or similar, and also write your own numeric solver instead of using ode15s( ). Probably not worth it. Do you really know those values to that precision anyway? Where do they come from?

请先登录，再进行评论。