Solving the two body problem in Matlab

[patric44]

Homework Statement

solve the two body problem using a Matlab script ?

Relevant Equations

a = -mu/r^3

hi guys
i am trying to code an algorithm for computing the trajectory of a basic two body problem situation according to the equation
$$\ddot{r} = \frac{-\mu}{r^3} \vec{r}$$
i am trying to use the Euler method , but the problem is in converting this problem into a 3 separate equations one for each component
of $r$ .
can i separate this equation for each component of $\ddot{r}$ as
$$\ddot{x} = \frac{-\mu}{r^3} \vec{x} \quad \ddot{y} = \frac{-\mu}{r^3} \vec{y} \quad \ddot{z} = \frac{-\mu}{r^3} \vec{z} \quad $$
and solve each part using euler method separately , or i cannot do that ?

 

[DrClaude]

can i separate this equation for each component of $\ddot{r}$ as
$$\ddot{x} = \frac{-\mu}{r^3} \vec{x} \quad \ddot{y} = \frac{-\mu}{r^3} \vec{y} \quad \ddot{z} = \frac{-\mu}{r^3} \vec{z} \quad $$
and solve each part using euler method separately , or i cannot do that ?

Yes you can, and should.

The usual approach is to reduce the system to a set of coupled first order differential equations, by introducing new variables corresponding to $\dot{x}$, etc. Note that the Euler method is unstable. You should try a Runge-Kutta method.

You can also directly integrate the second-order equations using Verlet integration.

 

[patric44]

Yes you can, and should.

The usual approach is to reduce the system to a set of coupled first order differential equations, by introducing new variables corresponding to $\dot{x}$, etc. Note that the Euler method is unstable. You should try a Runge-Kutta method.

You can also directly integrate the second-order equations using Verlet integration.

i tried to that but it only works in the first iteration from my initial vector i provide , i checked every thing but for some reason the velocity resulting from each solution individually is not correct , can some one revise my code :

clc
clear all
x0 = [6569.944337;-1.387381];   % inital condition  [x,v_x]
y0 = [1152.258896;6.658808];   % inital condition  [y,v_y]
z0 = [390.787865;3.658486];   % inital condition  [y,v_z]
r1 = [6569.944337;1152.258896;390.787865];
R = 6.371e6;    %%meters
M = 5.972e24;  %%kg
G = 6.67e-11;   %%% some SI unit
mu = G*M;
dt  = 60;   % time step
T =3*dt;    % amount of time to integrate
% the forward euler
xF(:,1) = x0; %  [.....] tow of them above each other
yF(:,1) = y0; 
zF(:,1) = z0;
rn(:,1) =[xF(1,1);yF(1,1);zF(1,1)]; %% all the elements in the first array contain r1
tF(1) = 0;
for k = 1:T/dt
     A = [0 1;-mu/(norm(rn(:,k)))^3 0];
     tF(k+1) = k*dt;    
     xF(:,k+1) = (eye(2)+A*dt)*xF(:,k);
     yF(:,k+1) = (eye(2)+A*dt)*yF(:,k);
     zF(:,k+1) = (eye(2)+A*dt)*zF(:,k);
rn(:,k+1) = [xF(1,k+1);yF(1,k+1);zF(1,k+1)];
end
rn

 

[Dr Transport]

that is a huge time step... and the total integration time is only 3 time steps. I'd start there. Euler is very unstable wrt time step size

 

[patric44]

that is a huge time step... and the total integration time is only 3 time steps. I'd start there. Euler is very unstable wrt time step size

i have no choice here i have an initial data set came from a simulation , the interval between the data is 60 sec
, i know it seems the main issue but i guess something else causes an error in the velocitiy components

 

[patric44]

The program seems that the velocity components calculated is way way off but the positions seems to be correct, so when the second loop starts it begins with the vectors containing the incorrect velocity and every thing goes wrong, can anyone help

 

[DrClaude]

i have no choice here i have an initial data set came from a simulation , the interval between the data is 60 sec

The integration time step does not have to correspond to the output time step.

The program seems that the velocity components calculated is way way off but the positions seems to be correct, so when the second loop starts it begins with the vectors containing the incorrect velocity and every thing goes wrong, can anyone help

As you are using Matlab, you should really convert the system to first-order equations and use ode45 to integrate.