Hello all.(adsbygoogle = window.adsbygoogle || []).push({});

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:

Where W is the function for the acceleration: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

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.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

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.

**Physics Forums | Science Articles, Homework Help, Discussion**

Dismiss Notice

Join Physics Forums Today!

The friendliest, high quality science and math community on the planet! Everyone who loves science is here!

The friendliest, high quality science and math community on the planet! Everyone who loves science is here!

# Error in N-Body simulation

**Physics Forums | Science Articles, Homework Help, Discussion**