ode45 code for 2 body problem

Hi,

I wrote a code to simulate the orbits of Sirus A & B.

As you can see, I get a Spiral and not an elipse.

What is the probelm?

Thank you

Main_TB()

function Main_TB

  global m1 m2

  m1=3.978e30; %sirius A mass in kg
  m2=2.025e30; %sirius B mass in kg

% about one year in seconds
  tspan = [0, 1.578e8];


% initial_condtions = [X1, Y1, U1, V1, X2, Y2, U2, V2]

X1 = -9.991e10;
Y1 = 0;
U1 = 0;
V1 = 11627.08;
X2 = 1.9629e11;
Y2 = 0;
U2 = 0;
V2 = -26231.69;

  initial_conditions = [-9.991e10, 0, 0, 11627.08, 1.9629e11, 0, 0, -26231.69];

  opt = odeset('AbsTol', 1e-8, 'RelTol', 1e-8);

  [t,z] = ode45('TB',tspan, initial_conditions);

  plot(z(:,1),z(:,3),"-");
  hold on;
  plot(z(:,5),z(:,7),"-");
  xlabel("X position (m)");
  ylabel("Y Position (m)");
  title ("Sirius A and B Orbits");
  axis equal; %makes sure circular elliptical orbits aren't distored
end  
  function dz = TB(t,z)
  global m1 m2
dz = zeros(8,1);
% Z(1)=X1
% Z(2)=U1 ( dx1/dt)
% Z(3)=Y1
% Z(4)=Y1 (dy1/dt)
% Z(5)=X2
% Z(6)=U2 (dx2/dt)
% Z(7)=Y2
% Z(8)= V2 (dy2/dt)

G=6.67e-11;


R12=sqrt((z(1) - z(5)).^2+(z(3) - z(7)).^2);

dz(1) = z(2);
dz(2) = (G*m2)*(z(5) - z(1))./R12.^3;
dz(3) = z(4);
dz(4) = (G*m2)*(z(7) - z(3))./R12.^3;
dz(5) = z(6);
dz(6) = (G*m1)*(z(1) - z(5))./R12.^3;
dz(7) = z(8);
dz(8) = (G*m1)*(z(3) - z(7))./R12.^3;
  end

[Torsten]

This is not a solution, but you forgot to include "opt" in the call to "ode45".

Answer 1

[James Tursa]

Your system has an initial non-zero total linear momentum. Recall that the 2-Body Problem has a constant center of mass motion integral as part of the solution. In fact, the N-Body Problem has this integral also. Your system is simply moving through inertial space about a constant linearly moving center of mass, hence the inertial plots look like spirals. There are various ways to accomodate this effect, such as changing to relative equations of motion instead of using inertial equations of motion as you are doing. That being said, I will simply subtract the center of mass from your solution to show the plot you probably expected to see with this constant inertial motion effect removed. I added a second plot just to show the linear motion of the center of mass as expected. And a third plot showing the velocity components of the center of mass. Ideally, the third plot should be exactly constant. The fact that it is not is just showing the numerical integration error.

Main_TB

function Main_TB
global m1 m2
m1=3.978e30; %sirius A mass in kg
m2=2.025e30; %sirius B mass in kg
% about one year in seconds
tspan = [0, 1.578e8];
% initial_condtions = [X1, U1, Y1, V1, X2, U2, Y2, V2] % CHANGED
X1 = -9.991e10;
Y1 = 0;
U1 = 0;
V1 = 11627.08;
X2 = 1.9629e11;
Y2 = 0;
U2 = 0;
V2 = -26231.69;
initial_conditions = [-9.991e10, 0, 0, 11627.08, 1.9629e11, 0, 0, -26231.69];
opt = odeset('AbsTol', 1e-8, 'RelTol', 1e-8);
[t,z] = ode45(@TB,tspan, initial_conditions,opt);
cm = (m1 * z(:,[1,3]) + m2 * z(:,[5,7])) / (m1 + m2); % Added center of mass calculation
z(:,[1,3]) = z(:,[1,3]) - cm; % Added
z(:,[5,7]) = z(:,[5,7]) - cm; % Added
figure % Added
plot(z(:,1),z(:,3),"-");
hold on;
plot(z(:,5),z(:,7),"-");
xlabel("X position (m)");
ylabel("Y Position (m)");
title ("Sirius A and B Orbits (Center of Mass Motion Removed)"); % CHANGED
axis equal; %makes sure circular elliptical orbits aren't distored
grid on % Added
% Added the following plots ...
figure
plot(cm(:,1),cm(:,2));
title("Sirius A and B Orbits - Center of Mass Motion")
xlabel("X position (m)");
ylabel("Y Position (m)");
axis equal;
grid on
figure
dcm = diff(cm) ./ diff(t);
tplot = t(1:end-1);
subplot(2,1,1);
plot(tplot,dcm(:,1));
title("Sirius A and B Orbits - Center of Mass Velocity")
xlabel("Time (s)");
ylabel("Xdot (m/s)");
grid on
subplot(2,1,2);
plot(tplot,dcm(:,2));
title("Sirius A and B Orbits - Center of Mass Velocity")
xlabel("Time (s)");
ylabel("Ydot (m/s)");
grid on
end  
function dz = TB(t,z)
global m1 m2
dz = zeros(8,1);
% Z(1)=X1
% Z(2)=U1 ( dx1/dt)
% Z(3)=Y1
% Z(4)=Y1 (dy1/dt)
% Z(5)=X2
% Z(6)=U2 (dx2/dt)
% Z(7)=Y2
% Z(8)= V2 (dy2/dt)
G=6.67e-11;
R12=sqrt((z(1) - z(5)).^2+(z(3) - z(7)).^2);
dz(1) = z(2);
dz(2) = (G*m2)*(z(5) - z(1))./R12.^3;
dz(3) = z(4);
dz(4) = (G*m2)*(z(7) - z(3))./R12.^3;
dz(5) = z(6);
dz(6) = (G*m1)*(z(1) - z(5))./R12.^3;
dz(7) = z(8);
dz(8) = (G*m1)*(z(3) - z(7))./R12.^3;
end

[Vered]

can you elborate what the equation should be? if we don't want to substract the center of mass?

[James Tursa]

E.g., depending on whether you want the motion relative to one of the bodies, or to the center of mass, you could look at these equations:

https://orbital-mechanics.space/the-n-body-problem/two-body-relative-motion.html

Answer 2

[Walter Roberson]

Use axis equal

Your x axis range is different from your y axis range, so you are seeing a distorted view.

[Vered]

Thank you, but my code has axis equal allready. so this is not the problem.

Can you think of something else ?