1. Limited time only! Sign up for a free 30min personal tutor trial with Chegg Tutors
    Dismiss Notice
Dismiss Notice
Join Physics Forums Today!
The friendliest, high quality science and math community on the planet! Everyone who loves science is here!

Error in N-Body simulation

  1. Jan 29, 2012 #1
    Hello all.

    I'm currently writing a program in Matlab that works out the positions of N bodies under the influence of gravity. The code is setup such that it only requires a mass and an initial position and velocity for each body. Any number of bodies can be entered but for testing purposes I am just using two.

    The code uses either Euler or RK4 (fourth-order Runge-Kutta integration). The problem I am getting occurs with either type of integration. Firstly, here is the code:

    Code (Text):
    clearvars

    % UNIVERSAL CONSTANTS

    G = 6.673*10^-11; % gravitational constant
    Msun = 1.98892*10^30;
    Mjup = 1.8987*10^27;
    Mearth = 5.9742*10^24;
    Mlunar = 7.349*10^22;
    AU = 1.49598*10^11;
    Vearth = 2.9786*10^4;
    Vlunar = 1.076*10^3;
    Tearth = 3.15576*10^7;
    Tmoon = 2.36058*10^6;

    % USER INPUTS

    IntegrationMethod = 2; % 1 = Euler, 2 = RK4
    N = 2; % The number of bodies
    m = [1*Msun 1*Mearth]; % Mass vector
    maxtime = 5*Tearth; % Maximum time to integrate over
    npts = 500; % Number of integration points

    % SETUP INTEGRATION

    H = maxtime/npts;

    % APPLY INITIAL CONDITIONS

    r_old(:,:,1) = [0 0 0]; % Sun starts at the origin
    r_old(:,:,2) = [0 AU 0]; % Earth starts 1 AU away

    drdt_old(:,:,1) = [0 0 0]; % Sun is initially motionless
    drdt_old(:,:,2) = [1*Vearth 0 0]; % Earth is moving perpendicular to the line joining it to the Sun

    % SETUP FIGURE

    close all
    figure(1)
    set(get(gca, 'xlabel'),'FontSize',20)
    set(get(gca, 'ylabel'),'FontSize',20)
    set(get(gca, 'title'),'FontSize',20)
    colourarray = ['r','b','g','g'];


    % DO INTEGRATION

    xlabel 'x'
    ylabel 'y'
    title 'N-Body Simulation'
    axis equal
    hold on


    for i = 1:npts
        if IntegrationMethod == 1
            drdt_new = drdt_old + H*W(N,G,m,r_old);
            r_new = r_old + H*drdt_new;
        elseif IntegrationMethod == 2
            k1 = zeros(1,3,N);
            k2 = zeros(1,3,N);
            k3 = zeros(1,3,N);
            k4 = zeros(1,3,N);
            l1 = zeros(1,3,N);
            l2 = zeros(1,3,N);
            l3 = zeros(1,3,N);
            l4 = zeros(1,3,N);
           
            k1 = drdt_old;
            l1 = W(N,G,m,r_old);
           
            k2 = drdt_old + H*l1/2;
            l2 = W(N,G,m,r_old + H*k1/2);
           
            k3 = drdt_old + H*l2/2;
            l3 = W(N,G,m,r_old + H*k2/2);
           
            k4 = drdt_old + H*l3;
            l4 = W(N,G,m,r_old + H*k3/2);
           
            drdt_new = drdt_old + (H/6)*(l1 + 2*l2 + 2*l3 + l4);
            r_new = r_old + (H/6)*(k1 + 2*k2 + 2*k3 + k4);
        end
        for j = 1:N
            scatter(r_new(1,1,j)/AU,r_new(1,2,j)/AU,colourarray(j),'.')
        end
        pause(0.0000001)
        r_old = r_new;
        drdt_old = drdt_new;
       
    end
    Where W is the function for the acceleration:

    Code (Text):
     function [accel] = W(N,G,m,r)
     accel = zeros(1,3,N);
     for j = 1:N
         for i = 1:N
             if i ~= j
                accel(:,:,j) = accel(:,:,j) + G*m(i)*(r(:,:,i)-r(:,:,j))/norm(r(:,:,i)-r(:,:,j))^3;
             end
         end
     end
     end
     
    Now here are the results. I have used Euler's method for these results but the RK4 results are pretty similar. RK4 is much more accurate than Euler integration but the error here is (I am pretty sure) not due to that. I have the Earth and Sun as the two bodies.

    http://img15.imageshack.us/img15/5015/2bodyorbitzoomedout.jpg [Broken]

    The above image looks pretty good yeah? The Earth orbits the Sun in pretty much a circular orbit. Now, let's zoom in on the Sun. Remember that the Sun also orbits the centre of mass of the Sun/Earth system but it will remain very close to the centre of mass because the Sun/Earth mass ratio is so large. I don't know about you though, but this doesn't look right:

    http://img171.imageshack.us/img171/3607/2bodyorbitzoomedin.jpg [Broken]

    I cannot explain this. Can someone please help? Just take the Euler integration...the code is so simple that I can't find the error. I have tried making the Sun body 1 and the Earth body 2, just in case I was referring to an index incorrectly or something but that did nothing.
     
    Last edited by a moderator: May 5, 2017
  2. jcsd
  3. Jan 29, 2012 #2

    edguy99

    User Avatar
    Gold Member

    I have not look closely at your code, but given the motion of the sun, I would guess what it would be. You are probably starting the earth to the right of the sun causing the sun to be attracted slightly to the right. When the earth is on the other side of the sun, it may well stop the motion of the sun to the right but thats all it does. You need to start the animation with the proper initial motion of the sun.

    Hope this helps...
     
  4. Jan 30, 2012 #3
    edguy99 is right, If you start the earth moving at speed v, and the sun motionless, there is an initial momentum of [itex] m_{earth} v [/itex], wich means you have a net drift of
    the center of mass by

    [tex] v \frac {m_{earth} } {m_{earth} + m_{sun}} [/tex]

    The orbit of the earth will drift as well, but the yearly drift is only

    [tex] 2 \pi r \frac { m_{earth} } {m_{earth} + m_{sun}} [/tex]

    about (1/50000) AU

    The graph produced for the sun is a cycloid, a superposition of a rotation round the center of mass and a linear motion with (nearly) the same speed. You get the same graph if you graph the position of a point on the circumference of a rolling wheel.
     
  5. Jan 30, 2012 #4
    Thank you kindly for your response. After analysing it a bit more I see exactly what you are saying. There isn't an error in the code, aside from the initial conditions which shouldn't be exactly zero for the initial velocity of the Sun.

    I have since tried some other examples too and have been able to replicate them (eg the figure 8 which is pretty cool).
     
Know someone interested in this topic? Share this thread via Reddit, Google+, Twitter, or Facebook