numerical simulation of Inverse law for orbital motion using rk4 method

Question

reema shrestha 2018-11-25

0
链接

此问题的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/431766-numerical-simulation-of-inverse-law-for-orbital-motion-using-rk4-method

编辑： madhan ravi 2018-11-25

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