Moving Object(s) Simulation

Hello all,

I started C++ programming a while back and quickly became unstuck after attempting my third, 'simple' program.

I am trying to simulate the path of three, interdependently moving masses through Cartesian space ((x,y) space if I'm not mistaken). Each point mass exerts a gravitational force on each of the others and is assigned with an initial velocity (speed and bearing). The desired result would be a text file containing the following information for each mass:
'(x,y) coordinates', 'speed', 'bearing' and 'elapsed time'

In a previous simulation, I was able to simulate one falling mass under the influence of another (which was fixed) by iterating the elapsed time value and reading off the previous iteration's values of speed and bearing. However, I was using fairly simple mathematics and when I tried to apply the same formulas to the new problem it felt like trying to smash a square block into a triangular hole.

Is anyone aware of the sort of mathematics needed for such a simulation? I've been pointed towards 'vector calculus' before but I'm not familiar with it - can anyone point me in the right direction?

Cheers!

How much (Maths) are you familiar with?
You can simply resolve the problem into x and y coordinates.... i.e. in terms of vectors .. i && j.
Then, the best and accurate way is to use calculus (just the basic matters..btw)...

Quote:

---

Originally Posted by manasij7479

How much (Maths) are you familiar with?
You can simply resolve the problem into x and y coordinates.... i.e. in terms of vectors .. i && j.
Then, the best and accurate way is to use calculus (just the basic matters..btw)...

---

I have done a fair bit of maths - A-level and now a little in Uni. This is mainly integration, optimization, differentiation, logic and a bit of matrix algebra though; I've never applied my mathematical knowledge to something like this.

I'm guessing that the use of vectors in this situation would involve matrix algebra? The last time I took a peek at vectors was roughly two years ago!

[P.s. Thanks for the quick reply :) ]

T
btw...Depending on your C/C++ knowledge , you could try making a maths library first..(I'm suggesting this because it is exactly what I'm doing now..)
before trying to simulate anything...

I'm not aware of how to create libraries. I don't suppose it's a file with pre-defined, mathematical functions...?

The basic idea is to compute the total force on each object due to all the other objects, divide the force by the object's mass to obtain an acceleration, then multiply the acceleration by the timestep to get a change in velocity, which you add to the current velocity, then multiply the velocity by the timestep to get a change in position. Do this over and over again.

If this isn't homework, and you want reference code, I can post some. To understand WHY it does what it does requires a bit of knowledge about physics, vectors, and differential equations.

Quote:

---

Originally Posted by brewbuck

The basic idea is to compute the total force on each object due to all the other objects, divide the force by the object's mass to obtain an acceleration, then multiply the acceleration by the timestep to get a change in velocity, which you add to the current velocity, then multiply the velocity by the timestep to get a change in position. Do this over and over again.

If this isn't homework, and you want reference code, I can post some. To understand WHY it does what it does requires a bit of knowledge about physics, vectors, and differential equations.

---

Cheers,
This project is something I want to work on over the summer (completely unrelated to my Economics course :D ), the above codes would be fantastic.

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.

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.

Quote:

---

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;
}

---

Quote:

---

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"

Fantastic!
That's plenty for my mind to chew on - this area of programming will keep me occupied for a fair while now :D

Cheers for the responses and the code, I'll post an application of the above when I get it all to work bug free.
Not too long into the distant future (if that makes any sense), I'll try the same problem but with bit-mapping!

Thanks again!

Out of curiosity, does this trade-off of between accuracy and data handling occur often in programming simulations?

In developing software, especially where geometry is involved, there's constantly a trade-off between accuracy and performance, despite that you're still always limited with respect to both and will spend a lot of time dealing with all kinds of boundary conditions.

One of the first things you may notice when you want to increase the accuracy of the simulation of an N-Body system, for example, has to do with collisions. Moving by discrete increments has this inevitable consequence of placing your object outside or inside of a boundary you mean to check for a collision against within a single frame's time. Typically that means within the single frame where you're moving from (t = time) x*t0 to x*t1 and say you're interested in the boundary which happens to be at x*t0 < x < x*t1, then you have to detect where along that path during that frame/iteration (x*t) the collision actually exists, from there you would have to move the object so it's in the 'proper' position at t1. I hope this example came through clearly, it's one of the more common kinds of problems you should expect to run in to.

Other problems with have to do with cramming the greatest number of particles that you can in to the simulation, interpolating/extrapolating over time for smoother visual motion, handling visual artifacts due to floating point precision limits, matrices vs quaternions for rotations, hardware acceleration, etc.

The discrete nature of computers introduces a lot of odd boundary conditions, and there is definitely a rich amount of knowledge to gain from those types of programming exercises, especially if you want to further develop your skill set in mathematics and logic.

The Euler method is the most simple implementation, it's also the least accurate and is prone to very quickly growing errors. As suggested you probably want to look at a higher order numerical integration, which mostly all work by interpolating each next iteration's movement by considering past and future points, things like that.

http://en.wikipedia.org/wiki/N-body_simulation
http://en.wikipedia.org/wiki/Euler_method
http://en.wikipedia.org/wiki/Runge%E...3Kutta_methods
http://en.wikipedia.org/wiki/Leapfrog_integration