ODE solver for restricted 3 body problem

조회 수: 8 (최근 30일)
Romio
Romio 2022년 9월 19일
답변: 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

채택된 답변

James Tursa
James Tursa 2022년 9월 19일
편집: James Tursa 2022년 9월 19일
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
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
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?

댓글을 달려면 로그인하십시오.

추가 답변 (1개)

Kevin Valencia
Kevin Valencia 2022년 12월 13일
This is the code working now, i am trying to do something similar, if it is correct i hope de plot be the answer that you needed!
Thanks for sharing!
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(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);
end

카테고리

Help CenterFile Exchange에서 Programming에 대해 자세히 알아보기

태그

제품


릴리스

R2021b

Community Treasure Hunt

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

Start Hunting!

Translated by