Dismiss Notice
Join Physics Forums Today!
The friendliest, high quality science and math community on the planet! Everyone who loves science is here!

Plotting Orbits on MATLAB

  1. Aug 16, 2016 #1
    I'm trying to plot the solutions of the second order differential equation d^2R/dt^2 = GM/R^2 + Lz^2/R^3. I'm reducing this to a system of first order ODEs and then using RK4 to solve this system.

    My code is given by
    Code (Text):
    function RK4system()                                      
    Tsim   = 10;           % simulate over Tsim seconds
    h      = 0.001;         % step size
    N      = Tsim/h;       % number of steps to simulate
    x4      = ones(2,N);
    t4      = zeros(1,N);
    theta  = ones(1,N);
    G      = 6.67*power(10,-11);
    M      = 1.5*1.989*power(10,41);
    Lz     = 2.623*power(10,2);
    R_1    = 2.623*power(10,20);
    U_1    = 0.5;
    x(1,1) = U_1;          %IC_ODE_1
    x(2,1) = (G*M/(power(R_1,2))+power(Lz,2)/(power(R_1,3)));  %IC_ODE_2

    x4 = ones(2,N);
    t4 = zeros(1,N);

    for n = 1:N-1
       t4(n+1)      = t4(n) + h;
       k1_4         = Orbit(t4(n), x4(:,n));
       k2_4         = Orbit(t4(n)+0.5*h, x4(:,n)+k1_4*0.5*h);
       k3_4         = Orbit(t4(n)+0.5*h, x4(:,n)+0.5*k2_4*h);
       k4_4         = Orbit(t4(n)+h, x4(:,n)+k3_4*h);
       phi4         = (1/6)*(k1_4 + 2*k2_4 + 2*k3_4 + k4_4);
       x4(:,n+1)    = x4(:,n) + phi4*h;
       theta(:,n+1)  = theta(:,n) + h*Angle(t(n), theta(:,n));
       xvar4(:,n+1) = x2(:,n+1)*cos(t4(n));
       yvar4(:,n+1) = x2(:,n+1)*sin(t4(n));

    plot(xvar4,yvar4, 'b-', 'LineWidth', 1)

    function dtheta = Angle(t,R)
    Lz = 100000*2.623*power(10,20);
    dtheta(1) = Lz/(power(R(1),2));

    function dx = Orbit(t,R)
    G           = 6.67*power(10,-11);
    M           = 1.5*1.989*power(10,41);
    Lz          = 100000*2.623*power(10,20);
    % Equations of motion
    dx          = zeros(2,1);
    dx(1)       = R(2);
    dx(2)       = (G*M/(power(R(1),2))+power(Lz,2)/(power(R(1),3)));
    I'm not getting the plot I want, If I set V = sqrt(GM/R) and Lz = R*sqrt(GM/R) I should get a circular orbit. If V > sqrt(GM/R) I should get a parabolic orbit. This is not what I'm getting.

    Please help
  2. jcsd
  3. Aug 16, 2016 #2
    Can you attach some plots for those cases?
    Images can help us spot some common errors more quickly than reading code.

    One obvious question, what does Orbit() mean?
    Finally I usually refrain from using functions like power and for example dot.
    Once you start running longer simulations this will slow down your program significantly due to checks built into those functions.
    In one case I was able to decrease the runtime of a simulation by a factor 10 by using ##a^\prime.a## for finding the length of a vector.
  4. Aug 16, 2016 #3

    Dr Transport

    User Avatar
    Science Advisor
    Gold Member

    did you try running this rk4 routine against something like a harmonic oscillator to check and see if it gives correct answers. I also suspect that the constants are giving it problems. R1 is on the order of [itex] 10^{20}[/itex] and when you square and cube it then divide, your essentially getting zero for some of the terms.
  5. Aug 16, 2016 #4

    Dr Transport

    User Avatar
    Science Advisor
    Gold Member

    [itex]T_{sim} [/itex] is only 10 seconds..... it takes a year to orbit the sun......., ya barely moved.
  6. Aug 16, 2016 #5


    User Avatar
    Education Advisor
    Gold Member

    Where are you defining t?
Share this great discussion with others via Reddit, Google+, Twitter, or Facebook

Have something to add?
Draft saved Draft deleted