numerical simulation of Inverse law for orbital motion using rk4 method

1 次查看(过去 30 天)
I a trying to simulate orbital motion of Earth and Sun using RK4 method. For r^2,I got an eliptical orbit but the velocities don't change as it should be according to kepler i.e. max at perihelion and min at apehelion,instead it has periodical change. I am not sure where I have got it wrong. How do I update the change in velocity?
Also, I tried change it to cube inverse law, for r^3, the orbit should have been unstable but it traces elliptical path.
Here is my code,
% t in years and distance in AU
%e = 0.15 tmax=1 B=2
function inverse_law_elliptical(e,tmax,B)
c=[0;1/2;1/2;1];
a=[0 0 0 0;1/2 0 0 0;0 1/2 0 0;0 0 1 0];
w=[1/6 1/3 1/3 1/6];
Ms=2*10^30; %mass of sun
G=4*pi^2/Ms; %6.67e-11*(6.68459e-12)^(3)*(3600*24*365)^2;
h=0.001; % step size
a0=(1/(0.998))^(1/3); %semi major axis of earth
xs=a0*e; %co-ordinates of sun
ys=0;
x0=a0; %initial position of earth
y0=0;
r0=sqrt((x0-xs)^2+y0^2); %initial distance between sun and earth
vx=0; %initial velocities of earth
vy=2*pi*r0;
t(1)=0; %initalizing t
i=1; %initializing i
s(:,1)=[x0;y0;vx;vy];
r=@(t) sqrt((s(1)-xs)^2+s(2)^2);
F=@(t,s) [s(3);s(4);-G*Ms*s(1)/r(t)^(B+1);-G*Ms*s(2)/r(t)^(B+1)]; %function
%initializing RK4 method
while t(i)<tmax
k=zeros(length(s(:,1)),length(c));
for j=1:length(c)
k(:,j)=h*F(t(i)+c(j)*h,s(:,i)+k*a(j,:)');
end
s(:,i+1)=s(:,i)+k*w';
t(i+1)=t(i)+h;
i=i+1;
end
x1=s(1,:); %x and y of earth
y1=s(2,:);
t=length(s(1,:));
% %animation loop
for i=1:t
plot(xs,ys,'.y','Markersize',70)
hold on
plot(x1(1:i),y1(1:i),'k','LineWidth',2)
hold off
axis([-2 2 -2 2]);
title('Simulation of Earth around Sun','fontsize',12)
xlabel('x-axis','fontsize',10)
ylabel('y-axis','fontsize',10)
pause(0.01)
end
  2 个评论
KALYAN ACHARJYA
KALYAN ACHARJYA 2018-11-25
编辑:madhan ravi 2018-11-25
You have to tell us which variable name represent the velocity, I can only see initial velocity
vx=0;
To change the velocity, you have to choose it as a vector or update in somewhete within the code.
Ask a specific question please!
reema shrestha
reema shrestha 2018-11-25
编辑:madhan ravi 2018-11-25
vx=0 is the initial velocity and s(3) is the position of velocity of x which is stored in s array.

请先登录,再进行评论。

回答(0 个)

类别

Help CenterFile Exchange 中查找有关 Earth and Planetary Science 的更多信息

标签

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by