Originally Posted by

**Distort**
I tried the above method shortly after my second project but found that the vector of the force acting on each mass changes with displacement; this might be solved for two freely moving masses (dealing with only one relative movement) but becomes inaccurate when working with three freely moving masses.

The problem lies in the continuous change in the vector of the force and the discrete nature of the time increment. Whilst this can be partly solved by making the time increments minutely small (equivalent to derivation - incurring a costly amount of data calculation), the calculations are still imperfect. Maybe I'm being too nitpickety and this sort of inaccuracy is required for any useful and feasible simulation - after all, it doesn't need to be inch-perfect, just representative.

You've recognized the core difficulty of this problem. There is no analytic solution for more than two bodies (this is the N-body problem). There are ways of numerically integrating the system which are more or less accurate. The naive method (leapfrog forward finite difference, the one I'll post here) requires a fairly small timestep to maintain accuracy. More complicated methods like Crank-Nicholson can use larger timesteps, but are difficult to implement. The time complexity is O(n^2) which makes it infeasible for simulating thousands of objects. In those cases a multipole approximation is used.

However, for the basic forward-difference method the calculations are not terribly intensive. Here's a skeleton you can start with (I do not guarantee it is bug free -- I whipped this up in the time it took you to reply to my post):

Code:

#include <cmath>
struct Body
{
double x, y;
double vx, vy;
double fx, fy;
double mass;
};
// Interact() calculates the force on a and b due to their mutual attraction
void Interact(double Gconst, // gravitational constant
Body *a, // first body to interact
Body *b) // second body to interact
{
// Compute displacement components in x and y directions
double dx = b->x - a->x;
double dy = b->y - a->y;
// Compute square of the distance, needed for force,
// and absolute distance, needed for vector projection
double r2 = dx * dx + dy * dy;
double r = std::sqrt(r2);
// Compute the force -- Newton's law of gravity.
//
// Check if r2 is zero (bodies are at the exact same location).
// This can't happen in the physical world, but it could happen
// in this simulation. If it does happen, let the force be zero.
// Note that if r2 is zero, then dx and dy are zero, thus the
// projection onto x and y directions will result in zero no matter
// what value is chosen for F in that case. Two bodies at the exact
// same point in space don't go exploding off in random directions,
// so this isn't unreasonable.
double F;
if (r2 != 0)
F = Gconst * a->mass * b->mass / r2;
else
F = 0;
// Project force onto dx and dy
double Fx = F * dx / r;
double Fy = F * dy / r;
// Update forces on particles -- forces are equal, opposite
a->fx += Fx;
a->fy += Fy;
b->fx -= Fx;
b->fy -= Fy;
}
// Update() computes the next position of the object for the given timestep
void Update(double dt, // timestep
Body *a) // body to update
{
// Update velocity
double scale = dt / a->mass;
a->vx += scale * a->fx;
a->vy += scale * a->fy;
// Update position
a->x += dt * a->vx;
a->y += dt * a->vy;
// Clear the force accumulators for the next iteration
a->fx = 0;
a->fy = 0;
}
void InteractAll(double Gconst, // gravitational constant
double dt, // timestep
Body *bodies, // array of bodies
int nBodies) // size of array
{
// For each pair of bodies...
for (int i = 0; i < nBodies - 1; ++i)
{
for (int j = i + 1; j < nBodies; ++j)
{
// Compute the force
Interact(Gconst, &bodies[i], &bodies[j]);
}
}
// TODO: Apply any external forces now
// ...
// For each body...
for (int i = 0; i < nBodies; ++i)
{
// Update the velocity and position
Update(dt, &bodies[i]);
}
}
int main()
{
// For demo, use three bodies
const int nBodies = 3;
Body bodies[nBodies];
// Make sure force accumulators start out at zero
for (int i = 0; i < nBodies; ++i)
{
bodies[i].fx = 0;
bodies[i].fy = 0;
}
// TODO: Set up the bodies' masses, initial positions, and initial velocities
//...
// Configure the gravitational constant and the timestep. Values are samples only,
// choosing the right value of dt is where some math knowledge comes in
const double Gconst = 1;
const double dt = 0.01;
// Simulate it
while (1)
{
InteractAll(Gconst, dt, bodies, nBodies);
// TODO: Print out whatever you want to print (position, velocity, etc)
}
return 0;
}

Ultimately, I want to use time increments to allow for the user to interject mid-simulation to act as an 'exogenous' shock (effectively 'knocking' one of the masses) and then see how it responds.

To implement such a thing, you'd apply an external force to the particles by tweaking the fx and fy at the spot in the code marked "Apply any external forces now"