Project 3: Orbital Mechanics

by Wyatt Rush, Sophaktra Suy, and John Wilson

Background

A good case study to examine the effect of solving differential equations is simulating the motion of an orbiting satellite. Keep in mind Newton's second law and the law of gravitation for two bodies of mass, giving us the equation $F = \frac{gm_1m_2}{r^2}$.

The One-Body Problem

As an introduction, we look at the one-body problem. In the one-body problem, one mass is considered negligible with the other, such as in the case of a small satellite orbiting a large planet. This offers us a simplified problem which allows us to neglect the force of the satellite on the planet so that the planet may be regarded as fixed.

To start, we place the large mass at the origin, and denote the position of the satellite by $(x,y)$. The distance between the mass and the satellite is $r = \sqrt{x^2+y^2}$ and the force on the satellite is in the direction of the large mass, so the direction vector is equal to $(-\frac{x}{\sqrt{x^2+y^2}}, -\frac{y}{\sqrt{x^2+y^2}})$. Taking all of this into account, the force on the satellite in terms of components is $(F_x, F_y) = (\frac{gm_1m_2}{x^2+y^2}\frac{-x}{\sqrt{x^2+y^2}}, \frac{gm_1m_2}{x^2+y^2}\frac{-y}{\sqrt{x^2+y^2}})$.

Inserting these forces into Newton's law of motion yields the second order equations $m_1x^{,,} = -\frac{gm_1m_2x}{(x^2+y^2)^{3/2}}$ and $m_1y^{,,} = -\frac{gm_1m_2x}{(x^2+y^2)^{3/2}}$. By introducing the variables $v_x = x^,$ and $v_y = y^,$, we can reduce these to a system of four first-order equations: $x^{,} = v_x, v_x^{,} = \frac{gm_2x}{(x^2+y^2)^{3/2}}, y^{,} = v_y, v_y^{,} = \frac{gm_2x}{(x^2+y^2)^{3/2}}$. While one might be inclined to use Euler's Method to approximate the solution to this, however, its practicality is severely limited in this case. Instead, we will replace the Euler step with the Trapezoid step for increased accuracy.

Although we use the one-body problem to lay the foundation for other problems, the one-body problem itself is essentially fictional since it ignores the force of the satellite on the larger mass. By factoring in the larger mass, we get what is called the two-body problem. We will address this one as well. After explaining that, we will discuss the three-body problem.

A Simpler Challenge: The Two-Body Problem

orbit.m

1234567891011121314151617181920212223242526272829303132333435363738394041424344454647484950515253545556575859%Program 6.4 Plotting program for one-body problem
% Inputs: time interval inter, initial conditions
% ic = [x0 vx0 y0 vy0], x position, x velocity, y pos, y vel,
% number of steps n, steps per point plotted p
% Calls a one-step method such as trapstep.m
% Example usage: orbit([0 100],[0 1 2 0],10000,5)
function z =orbit(inter,ic,n,p)
  v=VideoWriter('orbit10a','MPEG-4');open(v)
  h=(inter(2)-inter(1))/n; % plot n points
  x0=ic(1);vx0=ic(2);y0=ic(3);vy0=ic(4); % grab initial conds
  x1=ic(5);vx1=ic(6);y1=ic(7);vy1=ic(8); % ADDED A NEW BODY
  y(1,:)=[x0 vx0 y0 vy0 x1 vx1 y1 vy1]; t(1)=inter(1); % build y vector
  set(gca,'XLim',[-5 5],'YLim',[-5 5], 'XTick',[-5 0 5],'YTick',[-5 0 5]);
  head1=animatedline('color','r','Marker','.','markersize',35);
  tail1=animatedline('color','b','LineStyle','-');
  head2=animatedline('color','g','Marker','.','markersize',35);
  tail2=animatedline('color','c','LineStyle','-');
  %[px,py]=ginput(1); % include these three lines
  %[px1,py1]=ginput(1); % to enable mouse support
  %y(1,:)=[px px1-px py py1-py]; % 2 clicks set direction
  for k=1:n/p
     for i=1:p
         t(i+1)=t(i)+h;
         y(i+1,:)=trapstep(t(i),y(i,:),h);
     end
     y(1,:)=y(p+1,:);t(1)=t(p+1);
     clearpoints(head1);addpoints(head1,y(1,1),y(1,3))
     addpoints(tail1,y(1,1),y(1,3))
     clearpoints(head2);addpoints(head2,y(1,5),y(1,7))
     addpoints(tail2,y(1,5),y(1,7))
     drawnow;
     frame=getframe;
     writeVideo(v,frame);
  end
  close(v)

function y=eulerstep(t,x,h)
  %one step of the Euler method
  y=x+h*ydot(t,x);
function y = trapstep(t,x,h)
  %one step of the Trapezoid Method
  z1=ydot(t,x);
  g=x+h*z1;
  z2=ydot(t+h,g);
  y=x+h*(z1+z2)/2;
  function z = ydot(t,x)
  m1=0.03;m2=0.3;g=1;mg1=m1*g;mg2=m2*g;
  px1=x(1);py1=x(3);vx1=x(2);vy1=x(4);
  px2=x(5);py2=x(7);vx2=x(6);vy2=x(8);
  dist=sqrt((px2-px1)^2+(py2-py1)^2);
  z=zeros(1,8);
  z(1)=vx1;
  z(2)=(mg2*(px2-px1))/(dist^3);
  z(3)=vy1;
  z(4)=(mg2*(py2-py1))/(dist^3);
  z(5)=vx2;
  z(6)=(mg1*(px1-px2))/(dist^3);
  z(7)=vy2;
  z(8)=(mg1*(py1-py2))/(dist^3);

RK4.m

12345678function y = RK4(t,w,h)
    %one step of the Runge-Kutta order 4 method
    s1=ydot(t,w);
    s2=ydot(t+h/2,w+h*s1/2);
    s3=ydot(t+h/2,w+h*s2/2);
    s4=ydot(t+h,w+h*s3);
    y=w+h*(s1+2*s2+2*s3+s4)/6;
end

trapstep.m

1234567function y = trapstep(t,x,h)
    %one step of the Trapezoid Method
    z1=ydot(t,x);
    g=x+h*z1;
    z2=ydot(t+h,g);
    y=x+h*(z1+z2)/2;
end

ydot.m

12345678910111213141516171819202122232425262728293031323334353637function z = ydot(t,x)
    m1=.03;
    m2=.3;
    m3=.03;
    g=1;
    mg1=m1*g;
    mg2=m2*g;
    mg3=m3*g;
    px1=x(1);
    py1=x(3);
    vx1=x(2);
    vy1=x(4);
    px2 = x(5);
    py2 = x(7);
    vx2 = x(6);
    vy2 = x(8);
    px3 = x(9);
    py3 = x(11);
    vx3 = x(10);
    vy3 = x(12);
    dist12=sqrt((px2-px1)^2+(py2-py1)^2);
    dist13=sqrt((px3-px1)^2+(py3-py1)^2);
    dist23=sqrt((px3-px2)^2+(py3-py2)^2);
    z=zeros(1,12);
    z(1)=vx1;
    z(2)=(mg2*(px2-px1))/(dist12^3) + (mg3*(px3-px1))/(dist13^3);
    z(3)=vy1;
    z(4)=(mg2*(py2-py1))/(dist12^3) + (mg3*(py3-py1))/(dist13^3);
    z(5) = vx2;
    z(6)=(mg1*(px1-px2))/(dist12^3) + (mg3*(px3-px2))/(dist23^3);
    z(7)=vy2;
    z(8)=(mg1*(py1-py2))/(dist12^3) + (mg3*(py3-py2))/(dist23^3);
    z(9)=vx3;
    z(10)=(mg1*(px1-px3))/(dist13^3) + (mg2*(px2-px3))/(dist23^3);
    z(11)=vy3;
    z(12)=(mg1*(py1-py3))/(dist13^3) + (mg2*(py2-py3))/(dist23^3);
end

Outputs from orbit.m

Note: these files were generated using different output file names in the VideoWriter call on line 8 of orbit.m.

orbit9.mp4

Input: orbit([0 100], [2 0.2 2 -0.2 0 -0.02 0 0.02], 10000, 5), where $m_1 = 0.03, m_2 = 0.3$

orbit10a.mp4

Input: orbit([0 500], [0 0.1 1 -0.1 -2 -0.01 -1 0.01], 500000, 1000), where $m_1 = 0.03, m_2 = 0.3$

orbit10b.mp4

Input: orbit([0 500], [0 0.1 1 -0.1 -2 -0.01 -1 0.01], 500000, 1000), where $m_1 = 0.05, m_2 = 0.5$

orbit10c.mp4

Input: orbit([0 500], [0 0.1 1 -0.1 -2 -0.01 -1 0.01], 500000, 1000), where $m_1 = 0.08, m_2 = 0.8$

orbit11a.mp4

Input: orbit([0 500], [0 0.2 1 -0.2 -2 -0.2 -1 0.2], 500000, 1000), where $m_1 = m_2 = 2$

orbit11b.mp4

Input: orbit([0 500], [0 0.2 1 -0.2 -2 -0.2 -1 0.2], 500000, 1000), where $m_1 = m_2 = 1$

orbit11c.mp4

Input: orbit([0 500], [0 0.2 1 -0.2 -2 -0.2 -1 0.2], 500000, 1000), where $m_1 = m_2 = 0.5$

Addressing the Three-Body Problem

123456789101112131415161718192021222324252627282930313233343536373839404142function out = project_3(inter,ic,n,p)
    h = (inter(2) - inter(1))/n;
    x1 = ic(1);
    vx1 = ic(2);
    y1 = ic(3);
    vy1 = ic(4);
    x2 = ic(5);
    vx2 = ic(6);
    y2 = ic(7);
    vy2 = ic(8);
    x3 = ic(9);
    vx3 = ic(10);
    y3 = ic(11);
    vy3 = ic(12);
    y(1,:) = [x1 vx1 y1 vy1 x2 vx2 y2 vy2 x3 vx3 y3 vy3]; 
    t(1) = inter(1);
    set(gca,'XLim',[-5 5],'YLim',[-5 5],'XTick',[-5 0 5],'YTick',[-5 0 5]);
    head = animatedline('color','r','Marker','.','markersize',35);
    tail = animatedline('color','b','LineStyle','-');
    planet2 = animatedline('color','y','Marker','.','markersize',50);
    planet2tail = animatedline('color','g','LineStyle','-');
    planet3 = animatedline('color','b','Marker','.','markersize',35);
    planet3tail = animatedline('color','r','LineStyle','-');
    for k=1:n/p
        for i=1:p
            t(i+1)=t(i)+h;
            y(i+1,:)=RK4(t(i),y(i,:),h);
        end
        y(1,:)=y(p+1,:);
        t(1)=t(p+1);
        clearpoints(head);
        addpoints(head,y(1,1),y(1,3))
        addpoints(tail,y(1,1),y(1,3))
        clearpoints(planet2);
        addpoints(planet2,y(1,5),y(1,7))
        addpoints(planet2tail,y(1,5),y(1,7))
        clearpoints(planet3);
        addpoints(planet3,y(1,9),y(1,11))
        addpoints(planet3tail,y(1,9),y(1,11))
        drawnow;
    end
end

Outputs from project3.m

Note: these files were generated using different output file names in the VideoWriter call on line 8 of orbit.m.

q12.mp4

q13_10e-2.mp4

q13_10e-3.mp4

q13_10e-4.mp4

q16.mp4

q16_10e-1.mp4

q16_10e-3.mp4

References

Sauer, T. (2018). Numerical Analysis (3rd ed.) Boston: Pearson.