My attempt at an N-body simulation (gravity)

Over the previous year, I learnt about physics simulation and a lot about physics in general, mainly from brilliant. Thanks to those who helped me a lot along the way, Neeraj Anand Badgujar (Talulah Riley), Steven Chase, and Karan Chatrath.

First, we need data for each position of the $n$ bodies; you can do this by assigning random initial positions to the $n$ bodies within a certain range.

Build arrays that contain mass data of these bodies, $x$ and $y$ components of the position values, and $x$ and $y$ components of the velocities.

Now, here's how this brute-force CPU simulation works:

First, let's just try to figure out the dynamics of one body, for a small timestep.

Forces can be computed by multiplying the unit vectors to each other particle and the magnitudes of the gravitational force.

Essentially, 3 things should be computed each time for the $n-1$ bodies to find the total force exerted onto only 1 of these bodies.

The distance between the single body and each other body, the unit vectors, and finally gravitational force:

$\displaystyle \bold{F} = \frac{Gm_j m_i}{R^2} \hat{u}$

To do this as a program, loop through each of these bodies, and keep a total of the force which you keep updating until you reach the final $n-1$th body.

Then integrate using Euler integration.

This should be done to all the masses. Hence it will result in a nested loop. However, looping from 1 to the $n$th body results in a problem. When you start with the first body, you could be calculating the distance from the first body to the first body, which is zero. Then if you divide by the distance for the unit vector, this would result in division by zero.

So, you should only iterate the algorithm when the indexes of the position data $i$ and $j$ are not equal, ie. so you skip the the $j$th body when doing force computation.

Another consideration: what will happen if both bodies collide perfectly; ie. their centers touch? This would result in division by zero as well in the force computation; we cannot allow that to happen.

In my code I just let each body pass through one and other when they are a certain distance away from each other; this is done by setting the force to zero when they are a small distance away from each other, therefore they just move as per the previous values of their velocity.

However, once they are no longer in that small range, gravitational force again takes precedence.

If you want you can code in some resistive forces, like electrostatic force where all the charges are positive or all are negative, you can just calculate it along with the gravitational force using Coulomb's law. You can also add collision mechanics to it to prevent this error.

Also, when you want to animate each body, it will be a little computationally expensive. I found that using the Matlab rectangle() function to draw a filled circle less expensive than viscircles(). All you have to do is draw the rectangle with curvature and add a random colour to it. I personally added random colour using a colour matrix for random RGB values with row vectors for each body.

Like most of my simulations, I put the dynamics data into arrays which correspond with the time frame of the system and animate separately to the simulation.

However, since there are $n$ bodies, we should have 2 matrices $n \times f$ where $f$ is the number of frames the animation has, one for each component in $x$ and $y$.

The computation does run on the CPU, so wait a few seconds before the animation happens.

  1
  2
  3
  4
  5
  6
  7
  8
  9
 10
 11
 12
 13
 14
 15
 16
 17
 18
 19
 20
 21
 22
 23
 24
 25
 26
 27
 28
 29
 30
 31
 32
 33
 34
 35
 36
 37
 38
 39
 40
 41
 42
 43
 44
 45
 46
 47
 48
 49
 50
 51
 52
 53
 54
 55
 56
 57
 58
 59
 60
 61
 62
 63
 64
 65
 66
 67
 68
 69
 70
 71
 72
 73
 74
 75
 76
 77
 78
 79
 80
 81
 82
 83
 84
 85
 86
 87
 88
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122

clc;
clear;
close all;

time = 0;
deltaT = 10^-4;

G = 1;
n = 20;
totalForce = [0 0];

times = [];


%position data of bodies
xPositions = rand(1,n)*(10);
yPositions = rand(1,n)*(10);

%velocities
xVelocities = zeros(n,1);
yVelocities = zeros(n,1);

%Matrices for animation
xPositionMatrix = [];
yPositionMatrix = [];

%Mass data
Masses = [];

for i = 1:n
    Masses = [Masses; 5];
end

%% Simulation

count = 0;

while time <= 10
    if mod(count,100) == 0
        times = [times; time];
        xPositionMatrix = [xPositionMatrix; xPositions];
        yPositionMatrix = [yPositionMatrix; yPositions];
    end

    %% dynamics 

    %Outer loop (loop through each body, calculating dynamics of each body
    %before progressing the simulation
    for j = 1:n
        for i = 1:n
            if i~=j
                %Distance, unit vector, and total force (which is updated
                %every iteration
                currentDistance = hypot(xPositions(i)-xPositions(j),yPositions(i)-yPositions(j));
                unitVector = [(xPositions(i)-xPositions(j))/currentDistance (yPositions(i)-yPositions(j))/currentDistance];
                if currentDistance <= 0.3 %let bodies pass through each other
                    totalForce = [0 0];
                else %calculte gravitational force
                    totalForce = totalForce + unitVector*(G*Masses(j)*Masses(i)/currentDistance^2);
                end
            end
        end

        %% integrate and update positions for each body

        %Euler integration for jth body 
        xVelocities(j) = xVelocities(j) + totalForce(1)/Masses(j)*deltaT;
        yVelocities(j) = yVelocities(j) + totalForce(2)/Masses(j)*deltaT;

        xPositions(j) = xPositions(j) + xVelocities(j)*deltaT;
        yPositions(j) = yPositions(j) + yVelocities(j)*deltaT;

        %remember to reset total force to zero before redoing force
        %computation for a new jth body!
        totalForce = [0 0];
    end

    %progress the simulation and recalculate the dynamics for each body
    time = time + deltaT;
    count = count + 1;
end


%% animation
axis(gca, "equal");
axis([0 12 0 12]);
grid on;

particles = [];

particleColourMatrix = [];

r = 0.5;

for i = 1:n
    %random colour values of the particles
    particleColourMatrix = [particleColourMatrix; round(rand(1,3),1)];
end

for i = 1:length(times)
    title("time: " + num2str(times(i)) + "s");

    for j = 1:n
        %use rectangle() to draw a filled circle
        particles = [particles;     rectangle('Curvature', [1 1], ...
        'Position', [xPositionMatrix(i,j)-r yPositionMatrix(i,j) r r], ...
        'facecolor', particleColourMatrix(j,:), 'edgecolor', 'none') ];
    end

    pause(10^-6);


    if i < length(times)
        for j = 1:n
            %delete old elements before redrawing 
            delete(particles(j));
        end
    end

    %reset the array to zero before refilling the array
    particles = [];
end