# Gravitational Brute Force N-body Algorithm

I've just started dabbling with code, and as a learning exercise I've written a simple algorithm for solving gravitational n-body problems numerically in JavaScript.

I would be very grateful for general feedback as far as the quality of the code goes (first time I'm asking somebody to review my code; I'm very shy about it!), and more specifically advise on how I could optimize it as when I implement the algorithm in THREE.js in the shape of a simulation of the solar system the simulation starts lagging when I include more than around 70 bodies in the simulation. I know n-body algorithms are computationally intensive by nature, especially those of the brute force variety, but I'm pretty sure I've missed out on some useful tricks that would push it a bit further.

Update: to see a working example of the code, check out the jsfiddle I have thrown together https://jsfiddle.net/e9t88dwh

function nBodyProblem(parameters) {

this.g = parameters.g; //Gravitational constant
this.law = parameters.law; //Force law, inverse-square in our universe
this.dt = parameters.dt; //Time step to be used
this.masses = parameters.masses; //Object array with initial conditions for the bodiess to be simulated

}

nBodyProblem.prototype = {

constructor: nBodyProblem,

getDistance: function (m1, m2) {

return Math.pow((m1.x - m2.x) * (m1.x - m2.x) + (m1.y - m2.y) * (m1.y - m2.y) + (m1.z - m2.z) * (m1.z - m2.z), this.law);

},

updatePositionVectors: function () {

for (var i = 0, len = this.masses.length; i < len; i++) {

this.masses[i].x += this.masses[i].vx * this.dt;
this.masses[i].y += this.masses[i].vy * this.dt;
this.masses[i].z += this.masses[i].vz * this.dt;
}

return this;

},

updateVelocityVectors: function () {

var forceVectorX = 0;
var forceVectorY = 0;
var forceVectorZ = 0;

for (var i = 0, ilen = this.masses.length; i < ilen; i++) {

for (var j = 0, jlen = this.masses.length; j < jlen; j++) {

//We don't want to calculate self gravity!!!!!

if (i !== j) {

forceVectorX += (this.g * this.masses[j].m) * (this.masses[j].x - this.masses[i].x) / this.getDistance(this.masses[i], this.masses[j]);
forceVectorY += (this.g * this.masses[j].m) * (this.masses[j].y - this.masses[i].y) / this.getDistance(this.masses[i], this.masses[j]);
forceVectorZ += (this.g * this.masses[j].m) * (this.masses[j].z - this.masses[i].z) / this.getDistance(this.masses[i], this.masses[j]);

}

}

this.masses[i].vx += forceVectorX * this.dt;
this.masses[i].vy += forceVectorY * this.dt;
this.masses[i].vz += forceVectorZ * this.dt;

forceVectorX = 0;
forceVectorY = 0;
forceVectorZ = 0;

}

return this;

}

};


You didn't show the value you are using for parameters.law. Given that you used this in the function getDistance, I strongly suspect you are doing it wrong. You either shouldn't have called the function getDistance or you should have used a different equation. (The correct value for parameters.law given the expressions you are using is 1.5.)

What you should be using to compute the acceleration of a body due to gravitation is along the lines of

delta_x = this.masses[j].x - this.masses[i].x
delta_y = this.masses[j].y - this.masses[i].y
delta_z = this.masses[j].z - this.masses[i].z
distsq = delta_x*delta_x + delta_y*delta_y + delta_z* delta_z
dist_cubed = distsq * Math.sqrt(distsq)     # Note 1: Specialized to 3D space.
# Note 2: No call to pow.
fact = this.masses[j].mu / dist_cubed       # Note 3: Only one division.
# Note 4: Using mu rather than G*M
acceleration_x += delta_x * fact            # Note 5: Using previous calculations
# Note 6: It's acceleration, not force.
acceleration_y += delta_y * fact
acceleration_z += delta_z * fact


I'm going to go over those noted comments above in detail.

1. Specialized to 3D space.
Just as you can only tie your shoelaces in 3D space, gravity basically doesn't work in anything but 3+1 space (3 spatial dimensions, 1 time dimension).

2. No call to pow.
Unless you are using Fortran, it's best to avoid pow. Some compilers such as gcc and clang know how to deal with pow(x,2). Javascript? No chance. And even gcc and clang don't know how to deal with x3/2. Fortran compilers not only know how to efficiently compute x**2, x**3, they also know how to deal with things like x**(-2), x**(-3), and also x**(3/2).

3. Only one division.
In fact, you can do even better than what I did. There's no reason to compute the relation between particles j and i when you've already computed the relation be particles i and j.

4. Using mu (μ) rather than G*M.
We are taught in elementary physics that Newton's law of gravitation is F=GM1M2/R2. This is correct, but it is also a lie to children. The masses of the Sun, the planets, and anything else to which we humans have sent a probe are computed by dividing the observed gravitational parameter for that object by G. There's a problem here: the Newtonian gravitational constant G is by far the one physical constant that physicists cannot measure very accurately. It's best to avoid G, and there is no reason to use it. Use gravitational parameters instead. (Or treat mass as a derived unit, with base units of length3/time2. Same difference.)

5. Using previous calculations.
Your code is doing a lot of repeated calculations, and they are very expensive calculations. When you're doing scientific programming, it's important to pay attention to look out for redundant calculations.

6. It's acceleration, not force.
You made two key mistakes with regard to units. One was the calculation return Math.pow((m1.x - m2.x) * (m1.x - m2.x) + (m1.y - m2.y) * (m1.y - m2.y) + (m1.z - m2.z) * (m1.z - m2.z), this.law). You mislabeled this getDistance. This is not a distance function. The other is your calculation of force. Force has units of mass*acceleration. You are clearly using the quantity you calculated as acceleration. When you are doing scientific programming, it's extremely important to get your units right. While getting the units right does not guarantee that you have the physics right, getting the units wrong most certainly does guarantee that you have the physics wrong.

One last point: You did not mention how you are using updatePositionVectors and updateVelocityVectors. I'm assuming you used one of the following:

# Naive Euler integration.
updatePositionVectors()
updateVelocityVectors()


or

# Symplectic Euler integration.
updateVelocityVectors()
updatePositionVectors()


While neither one of these is a very good integrator, the latter choice is far better than the first. You should try the two variations and see what happens. Once you see that the second choice is far better, you might want to learn about more advanced integration techniques.

• Hi David! I have put together a jsfidlle where the inner solar system (sun, mercury, venus, earth and the moon) is simulated: jsfiddle.net/e9t88dwh - you can't see the moon because of the scale used, but it is there and appears to tag along with Earth in its orbit around the Sun without major issues. As far as I can tell, the physics seem to work. But yea, a leapfrog or rk iterator would be a step forward, but that would mean more computations, and yea.. Fewer bodies in the simulation! I guess my code is more of a toy than NASA's new algorithm for calculating spacecraft trajectories:) Feb 23 '16 at 13:15
• @HappyKoala - I forked your jsfiddle and made a few mods: jsfiddle.net/fuvf9xp8/3 . I switched to leapfrog. That's very easy to do. Advance velocity by half a step initially, and then regularly update position and velocity by full steps, in that order. This lets you take bigger steps. Change dt to 0.01 and the Moon stays coupled the Earth, but not in your simulation. I also switch to the scheme outlined in my answer. Finally, I offset the positions and velocities so as to make the center of mass of the system be at the origin and not be moving. Feb 23 '16 at 16:10
• Great stuff! I'm going to use your input in my 3D simulation of the solar system and I believe in giving credit when it is due, so would you like to contribute to the project on GitHub? If not I will just leave a reference to this discussion in the comments at the top of the file. One question, though: what do you mean by setting the center of mass to the origin of the grid by offsetting the positions and velocities? I thought the bary center of the solar system moves over time and that the bary center between the Sun and Jupiter is sometimes even above the surface of the sun? Feb 24 '16 at 10:55
• @HappyKoala -- In a heliocentric system, that's true. In a barycentric system, it's the Sun that moves and the barycenter that remains fixed. The latter is, ignoring the rest of the universe, an inertial frame. The former is not. You have the Sun moving, so you aren't using a heliocentric system. (That's a good thing.) I suspect you obtained your numbers from some site that gives the positions and velocities of the Sun and all of the planets in barycentric coordinates (and then perhaps converted to a system where the unit of length is the AU and the unit of time is one year). Feb 24 '16 at 12:11
• You left out the gas giants, and it's the gas giants that are responsible for almost all of the Sun's motion about the barycenter. By leaving those out, you have a system in which the barycenter is moving. There's nothing wrong with that, per se, but that means if you let your original version run for hours and hours, you would come back and find that the entire inner solar system has vanished. Feb 24 '16 at 12:13

you did something in updateVelocityVectors that I couldn't quite figure out at first

updateVelocityVectors: function () {

var forceVectorX = 0;
var forceVectorY = 0;
var forceVectorZ = 0;

for (var i = 0, ilen = this.masses.length; i < ilen; i++) {

for (var j = 0, jlen = this.masses.length; j < jlen; j++) {

//We don't want to calculate self gravity!!!!!

if (i !== j) {

forceVectorX += (this.g * this.masses[j].m) * (this.masses[j].x - this.masses[i].x) / this.getDistance(this.masses[i], this.masses[j]);
forceVectorY += (this.g * this.masses[j].m) * (this.masses[j].y - this.masses[i].y) / this.getDistance(this.masses[i], this.masses[j]);
forceVectorZ += (this.g * this.masses[j].m) * (this.masses[j].z - this.masses[i].z) / this.getDistance(this.masses[i], this.masses[j]);

}

}

this.masses[i].vx += forceVectorX * this.dt;
this.masses[i].vy += forceVectorY * this.dt;
this.masses[i].vz += forceVectorZ * this.dt;

forceVectorX = 0;
forceVectorY = 0;
forceVectorZ = 0;

}

return this;

}


in both of the for loops you added an extra variable to hold the length of this.masses, but never use it anywhere else in the loop. I think that you shouldn't declare the variable at all and just use this.masses.length in the declaration of both for loops like this

for (var j = 0, j < this.masses.length; j++) {


or an even better idea would be to assign the value to a variable outside of both loops so that you only call the length property of this.masses once, like this.

var massesLength = this.masses.length;
for (var i = 0; i < massesLength; i++) {
for (var j = 0; j < massesLength; j++) {
...


Same thing here as well

updatePositionVectors: function () {

for (var i = 0, len = this.masses.length; i < len; i++) {

this.masses[i].x += this.masses[i].vx * this.dt;
this.masses[i].y += this.masses[i].vy * this.dt;
this.masses[i].z += this.masses[i].vz * this.dt;
}

return this;

},


and I would remove some of the vertical white space as well.

updatePositionVectors: function () {
var massesLength = this.masses.length;
for (var i = 0; i < massesLength; i++) {
this.masses[i].x += this.masses[i].vx * this.dt;
this.masses[i].y += this.masses[i].vy * this.dt;
this.masses[i].z += this.masses[i].vz * this.dt;
}
return this;
},


@Ismael Miguel is correct, by using the length property inside the for declaration you would access that property every loop, which would make this loop very inefficient, so create the variable outside of the for loop.

• Thanks Malachi for your input! The idea behind the extra variable was indeed, as @IsmaelMiguel pointed out, to improve the performance of the loop by avoiding having to access the length property of this.masses every iteration. The code does look better when you declare the variable before the foor loop, though :)! Feb 23 '16 at 11:25