Brute Force N Body Implementation in C++

Question

I have wrote the following code in C++ for the n-body problem. This code is sequential as later on I am planning to parallelize it using OpenMP. I want to know whether I have wrote the correct implementation, if there are some errors or bugs, or if this is an invalid approach.

#include <iostream>
#include <ctime>
#include <math.h>
using namespace std;
#define N 100
#define G 6.673e-11
#define TIMESTAMP 1e11
struct Particle{
    double rx, ry;//position components
    double vx, vy;//velocity components
    double fx, fy;//force components
    double mass;//mass of the particle

};
Particle Update(Particle p, double timestamp)
{
    p.vx += timestamp*p.fx / p.mass;
    p.vy += timestamp*p.fy / p.mass;
    p.rx += timestamp*p.vx;
    p.ry += timestamp*p.vy;
    return p;
}
void PrintParticle(Particle p)
{
    printf("rx == %f ry == %f vx == %f vy == %f mass == %f\n", p.rx,p.ry,p.vx,p.vy,p.mass);
}
//Reset the forces on particle
Particle ResetForce(Particle p)
{
    p.fx = 0.0;
    p.fy = 0.0;
    return p;
}
//Add force to particle a by particle b
Particle AddForce(Particle a,Particle b)
{
    double EPS = 3E4;      // softening parameter (just to avoid infinities)
    double dx = b.rx - a.rx;
    double dy = b.ry - a.ry;
    double dist = sqrt(dx*dx + dy*dy);
    double F = (G * a.mass * b.mass) / (dist*dist + EPS*EPS);
    a.fx += F * dx / dist;
    a.fy += F * dy / dist;
    return a;

}

int main()
{
    Particle particles[N];
    srand(time(NULL));
    //randomly generating N Particles
    for (int i = 0; i < N; i++){
        double rx = 1e18*exp(-1.8)*(.5 - rand());
        particles[i].rx = rx;
        double ry = 1e18*exp(-1.8)*(.5 - rand());
        particles[i].ry = ry;
        double vx = 1e18*exp(-1.8)*(.5 - rand());
        particles[i].vx = vx;
        double vy = 1e18*exp(-1.8)*(.5 - rand());
        particles[i].vy = vy;
        double mass = 1.98892e30*rand()*10 + 1e20;
        particles[i].mass = mass;

    }

    int numberofiterations = 10;
    int count = 0;
    while (count < numberofiterations){
        for (int i = 0; i < N; i++)
        {
            particles[i] = ResetForce(particles[i]);
            for (int j = 0; j < N; j++)
            {
                if (i != j)
                {
                    particles[i] = AddForce(particles[i], particles[j]);
                }

            }
        }
        //loop again to update the time stamp here
        for (int i = 0; i < N; i++)
        {
            particles[i] = Update(particles[i], TIMESTAMP);
        }
        for (int i = 0; i < N; i++)
        {
            PrintParticle(particles[i]);
        }
        count++;
    }




    return 0;
}

Answer 1

Here are some observations that may help you improve your code.

Make sure to #include all required headers

This program calls printf and srand but does not include the corresponding headers. Fix that by adding these lines:

#include <cstdio>
#include <cstdlib>

Use objects

You have a Particle structure and then separate functions that operate on Particle data. With only a slight syntax change, you would have a real object instead of C-style code written in C++. For example, you have a loop right now that says this:

for (int i = 0; i < N; i++)
{
    particles[i] = Update(particles[i], TIMESTAMP);
}

With objects, and with C++11 or newer, you could instead write this:

for (auto &p : particles)
    p.update(TIMESTAMP);

All that would be required would be to make Update a member function.

Prefer stream I/O to printf

The printf function was a capable workhorse for many years, but the C++ iostream library is better in a number of regards. Although it's often more typing for the programmer initially, it's better because it has better type checking, less possibility for runtime overhead, and it fits well with the rest of C++. So I would probably change your PrintParticle to this:

friend ostream& operator<<(ostream& out, const Particle &p) {
    return out << 
        "rx == " << p.rx_ <<
        " ry == " << p.ry_ <<
        " vx == " << p.vx_ <<
        " vy == " << p.vy_ <<
        " mass == " << p.mass_ << '\n';
}

Use const where practical

There are a number of things in your code that look and act like constants, such as numberofiterations in main(). Declare them as const for a program more resistant to error. Also, move the constants, such as N into main rather than having them at file scope.

Create a function rather than repeating code

The functions which create the particle use the same rand()-based function several times. Instead of repeating it, make it into a function:

double xrand() {
    return 1e18*exp(-1.8)*(.5 - rand());
}

Don't abuse using namespace std

Putting using namespace std at the top of every program is a bad habit that you'd do well to avoid. It isn't necessarily wrong to use, but be aware of when you absolutely shouldn't do it (such as in header files).

Consider using a better random number generator

If you are using a compiler that supports at least C++11, consider using a better random number generator. In particular, instead of rand, you might want to look at std::uniform_real_distribution and friends in the <random> header.

Consider altering your algorithm

At the moment, the force on each particle from every other particle is calculated exhaustively, but Newton tells us that the force actually acts equally (but opposite) on both particles, so your loop could be simplified and calculations speeded if you take advantage of that fact to only calculate force vectors once and then apply them to both.

Double check your constant values

It seems to me that 3e4 is a somewhat large value for EPS. Additionally, in the calculation of the random mass, the math is currently 1.98892e30*rand()*10 + 1e20, but why not simply write 1.98892e31*rand() + 1e20 instead?

Eliminate return 0 at the end of main

When a C++ program reaches the end of main the compiler will automatically generate code to return 0, so there is no reason to put return 0; explicitly at the end of main.

Answer 2

All points by @Edvard are good. Here are a couple more.

Encapsulation

All your functions mutate a Particle structure. So why not just make this methods on your class?

struct Particle{
    double rx, ry;//position components
    double vx, vy;//velocity components
    double fx, fy;//force components
    double mass;//mass of the particle

};
Particle Update(Particle p, double timestamp);
void PrintParticle(Particle p);
Particle ResetForce(Particle p);
Particle AddForce(Particle a,Particle b);

I would have written this as:

classe Particle
{
    // Members are all private.
    // You mutate and access them via methods.
    double rx, ry;//position components
    double vx, vy;//velocity components
    double fx, fy;//force components
    double mass;//mass of the particle

    public:

        Particle(double rx, double ry, double vx, double vy, double fx, double fy, double mass)
            :rx(rx) ,ry(ry)
            ,vx(vx) ,vy(vy)
            ,fx(fx) ,fy(fy)
            ,mass(mass)
        {} 

        // If these methods return `void` or `Particle&` depends
        // on context and usage patterns. I have left then as close
        // to the original functions as possible.
        Particle& update(double timestamp);
        void printParticle(std::ostream& str);
        Particle& resetForce();
        Particle& addForce(Particle b);

        friend std::ostream& operator<<(std::ostream& str, Particle const& data)
        {
            printParticle(str);
            return str;
        }
};

Don't use macros

#define N 100
#define G 6.673e-11
#define TIMESTAMP 1e11

Macros don't obey scope rules and have no type information (they are simple text substitutions). Prefer to use const variables (or constexpr in C++11);

constexpr int    N         = 100;
constexpr double G         = 6.673e-11;
constexpr double TIMESTAMP = 1e11;

Prefer Standard containers

Arrays are limited. They don't work well with a lot of other constructs in the language as they quickly degrade into pointers just by looking at them funnily (them people start worrying about memory management).

Particle particles[N];

It is better to use modern constructs like std::vector<> or std::array<> as these object maintain there type information correctly.

std::array<Particle, N>  particles;   // fixed size array of particles.
                                      // useful if you they can be 
                                      // build by default contructor

Alternatively:

 std::vector<Particle>   particles;   // Dynamically sized array.
 particles.reserve(N);                // Useful when built dynamically
                                      // from a loop as you can call
                                      // the constructor specifically.

Prefer for(;;) over while() when counter is used

int count = 0;
while (count < numberofiterations){
  // STUFF
  count++;
}

// Easier to read and reason (as the increment is obvious)
for(int count = 0; count < numberofiterations; ++count) {
    // STUFF
}

Standard algorithms

Look at the standard algorithms and use them when you can. Alternatively use the for() syntax for looping over containers.

Refactoring the main now looks like this:

int main()
{
    srand(time(NULL));

    std::vector<Particle> particles;
    particles.reserve(N);

    //randomly generating N Particles
    for (int i = 0; i < N; i++)
    {
        particles.emplace_back(getRand(), getRand(), 
                               getRand(), getRand(),
                               getRand(), getRand(),
                               getRandMass());
    }

    int constexpr numberofiterations = 10;

    for(int count = 0; count < numberofiterations; ++count) {
        for (int i = 0; i < N; ++i) {
            particles[i].resetForce();
            for (int j = 0; j < N; ++j) {
                if (i != j) {
                    particles[i].addForce(particles[j]);
                }
            }
        }
        //loop again to update the time stamp here
        for (auto& particle: particles) {
            particle.update(TIMESTAMP);
        }
        std::copy(std::begin(particles), std::end(particles),
                  std::istream_iterator<Particle>(std::cout));
    }
}

Answer 3

All the advice from @Edvard and @Loki Astari are good. However, one thing I'd also recommend you on doing is to use dimensionless units. That way you have two advantages:

• It is easier to generalize the system;
• You are less prone to loss of precision.

@Snowbody (fairly) asked if such system as for \$G=1\$. Yes, but not only. One should set the whole system of units.

Say, if you have the following energy:

\$ U = G \sum_{ij} \frac{m_i m_j}{r_{ij}}\$

You can define dimensionless units as

\$ m' = m/m_0 \$;

\$ r' = r/r_0 \$;

\$ U' = U/E_0 \$.

Thus,

\$U' = \frac{U}{E_0} = \frac{G}{E0} \sum_{ij} \frac{m_i m_j}{r_{ij}} \$.

Replacing \$ r\$ and \$ m \$ in \$ U' \$.

\$U' = \frac{U}{E_0} = \frac{G}{E_0} \sum_{ij} \frac{m_i' m_j' m_0^2}{r_{ij'} r_0} \$.

\$U' = \frac{G m_0^2}{E_0r_0} \sum_{ij} \frac{m_i' m_j'}{r_{ij}'} \$.

Now comes the cool part. We want the first term to unity. So, we say that

\$ \frac{G m_0^2}{E_0r_0} =1 \$

or

\$ \frac{G m_0^2}{r_0} = E_0 \$

\$U' =\sum_{ij} \frac{m_i' m_j'}{r_{ij}} \$.

Now one knows how to link the real values with the simulations. Say one wants to simulate the earth-moon. Then \$ m_0 \$ can be the mass of the Earth. Then, in the program, \$ m_{Earth} = 1.0\$ and \$ m_{Moon} = 0.012300 \$*. One can do the same thing for \$ r_0 \$.

* Moon mass in Earth units according to wikipedia

Answer 4

• why are you storing the force component instead of the acceleration? Do you expect particles' masses to change? Division is usually more expensive than multiplication, so do it as little as possible.
• Why are you calculating distance squared, then square rooting it, then squaring it again? The compiler may be able to optimize this one though.