Prob with MATLAB code (Planetary motion)

[cloudy14]

Hi,

I am writing a code for planetary motion using ODE with runge-kutta.

Currently it has only two body, however, I could not get the initial values to fit in properly, the bodies will not move in circular orbital motion. They are moving towards each other.
at the v=[10 0;-10 0], whenever I try to add a value for the vy, the second body will disappear completely from the screen. Help please!

function orbit
global X v n G m       % Make these params. global.

n=2;    %No. of planets
G = 4*pi^2; % Grav. const
m=[4 4];  %m=[400 1.23e-2]; 
X=[0 0;50 0]; %initial positions
v=[10 0;-10 0]; %initail velocities %17.9979
state=zeros(n,4); %Empty the state memory
for k=1:n, %Forming the state matrix from the previous values
    state(k,:)=[X(k,1) X(k,2) v(k,1) v(k,2)];
end

nStep = 2000; %number of steps
tau = 0.01; %time step (yr)
time = 0;

for iStep=1:nStep
    for i=1:n
        state(i,:) = rk4(state(i,:),time,tau,@dEqs,G,i);
        %update new values of X and v
        X(i,:)=[state(i,1) state(i,2)];
        v(i,:)=[state(i,3) state(i,4)];
    end
time = time + tau;
xpos1(iStep)=X(1,1);
ypos1(iStep)=X(1,2);
xpos2(iStep)=X(2,1);
ypos2(iStep)=X(2,2);
%xpos3(iStep)=X(3,1);
%ypos3(iStep)=X(3,2);
clf('reset');
hold on;
axis on
axis equal
axis([-150 150 -150 150]);
plot(X(1,1),X(1,2),'bo','MarkerFaceColor','r','MarkerSize',12);
plot(X(2,1),X(2,2),'bo','MarkerFaceColor','y','MarkerSize',12);
%plot(X(3,1),X(3,2),'bo','MarkerFaceColor','r','MarkerSize',12);
plot(xpos1(1:iStep),ypos1(1:iStep));
plot(xpos2(1:iStep),ypos2(1:iStep));
%plot(xpos3(1:iStep),ypos3(1:iStep));

pause(0.001)


%  clf reset
%    set(gcf,'numbertitle','off','name',' n-body')
%    n = size(state,1);
%    s = max(sqrt(diag(state*state')));
%    if n <= 3, s = 2*s; end
%    axis([-s s -s s])
%    axis square

end


function F = dEqs(x,t,G,num) % First-order differential
global n m X

a = zeros(1,2); 
tiny=1e-20;
    for j=1:n,
        if j~=num;
        a(1) =G* (a(1)+ m(j)* (X(j,1)-x(1,1))/norm((X(j,1)-x(1,1))+tiny)^3 );%x component
        a(2) =G* (a(2)+ m(j)* (X(j,2)-x(1,2))/norm((X(j,2)-x(1,2))+tiny)^3 );%x component        
        end
    end

F = [x(1,3) x(1,4) a(1) a(2)]; 

% Runge-Kutta integrator (4th order)
function xout = rk4(x,t,tau,derivsRK,param,count)
half_tau = 0.5*tau;
F1 = feval(derivsRK,x,t,param,count);
t_half = t + half_tau;
xtemp = x + half_tau*F1;
F2 = feval(derivsRK,xtemp,t_half,param,count);
xtemp = x + half_tau*F2;
F3 = feval(derivsRK,xtemp,t_half,param,count);
t_full = t + tau;
xtemp = x + tau*F3;
F4 = feval(derivsRK,xtemp,t_full,param,count);
xout = x + tau/6.*(F1 + F4 + 2.*(F2+F3));
return;

 

[Valkarie]

It looks like you are trying to use the Runge-Kutta method to simulate the motion of two planets around a central point. When you set the initial velocities of the two planets to [10 0; -10 0], they will be moving towards each other and it is not possible to get them to move in a circular orbit. In order for two planets to move in a circular orbit around each other, their initial velocities must be such that the centripetal force balances the gravitational force. This can be calculated using the equations of motion. The acceleration of the two planets due to gravitational forces can be calculated by: a_1 = F_grav(m_2 / r^2)a_2 = F_grav(m_1 / r^2)Where F_grav is the gravitational force between the two planets and r is the distance between them. The centripetal acceleration of the two planets is given by: a_1 = v_1^2 / ra_2 = v_2^2 / rWhere v_1 and v_2 are the velocities of the two planets. The initial velocities of the two planets can then be calculated by setting the gravitational and centripetal accelerations equal to each other and solving for the velocities: v_1 = sqrt(G*m_2 / r)v_2 = sqrt(G*m_1 / r)Once you have the initial velocities, you can use the Runge-Kutta method to simulate the orbital motion of the two planets. I hope this helps.