Using Matlab ode45 to solve a simple 2 body problem

[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 :

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

 

[pasmith]

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.

 

[patric44]

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 :

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

 

[pasmith]

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

For example, to solve the Lorenz system I would use

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:

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 <br /> \mathbf{R} = \frac{\mathbf{r}}{\|\mathbf{r}(0)\|}, \qquad<br /> T = t\left( \frac{\mu}{\|\mathbf{r}(0)\|^3} \right)^{\frac 12} which has the advantage that \mu = 1: <br /> \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]

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

For example, to solve the Lorenz system I would use

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:

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 <br /> \mathbf{R} = \frac{\mathbf{r}}{\|\mathbf{r}(0)\|}, \qquad<br /> T = t\left( \frac{\mu}{\|\mathbf{r}(0)\|^3} \right)^{\frac 12} which has the advantage that \mu = 1: <br /> \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 ?!

 

[pasmith]

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.