# Using Matlab ode45 to solve a simple 2 body problem

• patric44
In summary, the conversation discusses a problem with solving a two-body problem in MATLAB using a table of data from a simulation. The problem involves using the two body equation and separating it into three component equations, but previous methods such as Euler method and ode45 have not been successful. The solution may involve solving a single six-dimensional problem and defining a function to evaluate dy/dt.

#### patric44

Homework Statement
trying to solve the two body problem using matlab ode45
Relevant Equations
\ddot{x} = -mu*x/r^3
hi guys
i am having some truble solving a simple 2body problem in MATLAB , basicaly i have a table of data that came from a simulation of a satellite orbiting the Earth with 1 min interval between them , i will only take the first inital position and velocity from that table and i sould construct a MATLAB script the will predict the future positions and velocities .
the problem relied basicaly on the two body equation :
$$\ddot{r} = \frac{-\mu}{r^{3}}*\vec{r}$$
i separated this equation into three different component equations , last time i tried to use Euler method but some one pointed that the step size might be very large for Euler the euler method , it was working realivialy accurate for getting the position component but the velocity component is way off by a factor of ##10^{5}## which is very srange ?! , so i tried to use ode45 but with no success :
Matlab:
R = 6.371e6;    %%meters
M = 5.972e24;  %%kg
G = 6.67e-11;   %%% some SI unit
mu = G*M;
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]
r0 = [6569.944337;1152.258896;390.787865];
rn(:,1) = r0;
dt = 60;
T = 10*dt;
% Integration:
for k = 1:T/dt
A = [0 1;-mu/(norm(rn(:,k))^3) 0];
[t, X]  = ode45(@(t,x) A*x, 0:dt:T,x0);
[t, Y]  = ode45(@(t,y) A*y, 0:dt:T,y0);
[t, Z]  = ode45(@(t,z) A*z, 0:dt:T,z0);
rn(:,k+1) = [X(k+1,1);Y(k+1,1);Z(k+1,1)];
end
rn
any help will be appreciated

Last edited by a moderator:
I think you need to solve a single six-dimensional problem for $(x,y,z,v_x,v_y,v_z)$, not three separate two-dimensional problems. ode45 will perform multiple time steps to get the final value, rather than only one, so if you're only updating $x$ then of course your value for $r^2$ will get more and more inaccurate.

The problem is not linear; there is no need to calculate $dx/dt$ using a matrix multip[lication. Instead define a function which takes time and the six-dimensional state vector as arguments and returns the six-dimensional vector $(v_x, v_y, v_z, \dot v_x, \dot v_y, \dot v_z)$.

I would suggest looking at the van der Pol example on the ode45 reference page and adapting that to your case.

pasmith said:
I think you need to solve a single six-dimensional problem for $(x,y,z,v_x,v_y,v_z)$, not three separate two-dimensional problems. ode45 will perform multiple time steps to get the final value, rather than only one, so if you're only updating $x$ then of course your value for $r^2$ will get more and more inaccurate.

The problem is not linear; there is no need to calculate $dx/dt$ using a matrix multip[lication. Instead define a function which takes time and the six-dimensional state vector as arguments and returns the six-dimensional vector $(v_x, v_y, v_z, \dot v_x, \dot v_y, \dot v_z)$.

I would suggest looking at the van der Pol example on the ode45 reference page and adapting that to your case.
i tried to use a 6 *1 vector but also no success :
Matlab:
y = [6569.944337;-1.387381;1152.258896;6.658808;390.787865;3.658486];
R = 6.371e6;    %%meters
M = 5.972e24;  %%kg
G = 6.67e-11;   %%% some SI unit
mu = G*M;
r0 = [6569.944337;1152.258896;390.787865];
rn(:,1) = r0;
dt = 60;
T = 5*dt;
% Integration:
for k = 1:T/dt
dy = [y(2);-mu*y(1)/(norm(rn(:,k)))^3;y(4);-mu*y(3)/(norm(rn(:,k)))^3;y(6);-mu*y(5)/(norm(rn(:,k)))^3];
[t, Y]  = ode45(dy,0:dt:T,y);
rn(:,k+1) = [Y(k+1,1);Y(k+1,3);Y(k+1,5)];
end

You really need to define a function to evaluate dy/dt.

For example, to solve the Lorenz system I would use

Matlab:
function dydt = lorenz(t, y, sigma, rho, beta)
dydt = zeroes(3,1);
dydt(1) = sigma * (y(2) - y(1));
dydt(2) = y(1) * (rho - y(3)) - y(1);
dydt(3) = y(1) * y(2) - beta * y(3);
end

And then to solve it, I just need one call to ode45:

Matlab:
sigma = 10;
beta = 8/3;
rho = 28;
dt = 0.01;
T = 40;
y0 = [1;1;1];
[t, y] = ode45( @(t, y) lorenz(t, y, sigma, rho, beta), 0:dt:T, y0);

I would also suggest rescaling your variables to something appropriate such as $$\mathbf{R} = \frac{\mathbf{r}}{\|\mathbf{r}(0)\|}, \qquad T = t\left( \frac{\mu}{\|\mathbf{r}(0)\|^3} \right)^{\frac 12}$$ which has the advantage that $\mu = 1$: $$\frac{d^2\mathbf{R}}{dT^2} = - \frac{ \mathbf{R} }{\|\mathbf{R}\|^3}.$$ It also means that you aren't simulating 60 seconds worth of a system whose orbital period is measured in days, months, or years.

• • patric44 and jim mcnamara
pasmith said:
You really need to define a function to evaluate dy/dt.

For example, to solve the Lorenz system I would use

Matlab:
function dydt = lorenz(t, y, sigma, rho, beta)
dydt = zeroes(3,1);
dydt(1) = sigma * (y(2) - y(1));
dydt(2) = y(1) * (rho - y(3)) - y(1);
dydt(3) = y(1) * y(2) - beta * y(3);
end

And then to solve it, I just need one call to ode45:

Matlab:
sigma = 10;
beta = 8/3;
rho = 28;
dt = 0.01;
T = 40;
y0 = [1;1;1];
[t, y] = ode45( @(t, y) lorenz(t, y, sigma, rho, beta), 0:dt:T, y0);

I would also suggest rescaling your variables to something appropriate such as $$\mathbf{R} = \frac{\mathbf{r}}{\|\mathbf{r}(0)\|}, \qquad T = t\left( \frac{\mu}{\|\mathbf{r}(0)\|^3} \right)^{\frac 12}$$ which has the advantage that $\mu = 1$: $$\frac{d^2\mathbf{R}}{dT^2} = - \frac{ \mathbf{R} }{\|\mathbf{R}\|^3}.$$ It also means that you aren't simulating 60 seconds worth of a system whose orbital period is measured in days, months, or years.
ok i will try that but i have a little question shouldn't i make a for loop the will take the values produced from the ode45 and update the values of the radius vector ##r## . or that's not the case ?!

patric44 said:
ok i will try that but i have a little question shouldn't i make a for loop the will take the values produced from the ode45 and update the values of the radius vector ##r## . or that's not the case ?!

It's not necessary; what you will get from ode45 is essentially (1) an N x 1 array of time values at wnich the solution was calculated and (2) an N x 6 array of the solution value at those times. (You actually get slightly more than that; the details are given on the reference page I linked in my first post.)

You may learn more by actually implementing the fourth-order Runge-Kutta method yourself; that will require a for loop and you will need to evaluate the right hand side of the ODE four times in each step at different values of $t$ and $y$.