- #1
patric44
- 308
- 40
- 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 :
any help will be appreciated
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
Last edited by a moderator: