# ODE solver for restricted 3 body problem

7 views (last 30 days)

Show older comments

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 Comments

### Accepted Answer

James Tursa
on 19 Sep 2022

Edited: James Tursa
on 19 Sep 2022

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 Comments

James Tursa
on 19 Sep 2022

Edited: James Tursa
on 19 Sep 2022

### More Answers (0)

### See Also

### Categories

### Community Treasure Hunt

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

Start Hunting!