Intro to Gravitational Simulations

[Sagekilla]

Hi all, I recently was posting in another thread on how to properly integrate the equations of motion in gravitational simulations. I thought this would be a good opportunity to provide some basic knowledge on simulations, numerical differential equations, and a little bit of programming.

I'm writing this as a resource for anyone who'd like to learn how to simulations and provide a crash course, so to speak, on setting up some basic simulations.

The tutorial in full is posted in the below three posts. Enjoy!Edit: D'oh, I just realized I posted this in the wrong forum. Can a moderator kindly move this to the "Programming" section?

 

[Sagekilla]

Finding the Gravitational Interaction

To start with, I'll be talking about gravity. As is known, gravity is an attractive force between two masses. It's proportional to each of their masses, and inversely proportional to the square of the distance between the two objects.

Gravity is also proportional to some proportionality constant, "G", known as the Gravitational Constant.

In the form of an equation, we have:
$$(1.1): F_{g} = \frac{-G M_{i} M_{j}}{r^{2}}$$

Equation 1.1 is easily extended to normal, 3D space through the use of vectors:
$$(1.2): \vec{F_{g}} = \frac{-G M_{i} M_{j}}{r^{2}} \hat{r_i}$$

Where the vector r in equation 1.2 is:
$$(1.3): \vec{r} = \vec{r_{i}} - \vec{r_{j}}$$

And the scalar r is the usual length or norm of the vector r:
$$(1.4): r = |\vec{r}| = \sqrt{r_{x}^{2} + r_{y}^{2} + r_{z}^{2}}$$

And finally, the unit vector r-hat is simply vector r normalized to length 1:
$$(1.5): \hat{r} = \frac{\vec{r}}{|\vec{r}|}$$

But how is this useful at all? Recall Newton's second law, which states that the force exerted on an object is equal to it's mass times the acceleration:
$$(1.6): \vec{F} = m \vec{a}$$

By equating equations 1.2 and 1.6, we find:
$$(1.7): \vec{F} = m \vec{a} = \frac{-G M_{i} M_{j}}{r^{2}} \hat{r}$$

If we take the left hand side to be force on object i, then we simply have:
$$(1.8): M_{i} \vec{a_{i}} = \frac{-G M_{i} M_{j}}{r^{2}} \hat{r}$$
$$(1.9): \vec{a_{i}} = \frac{-G M_{j}}{r^{2}} \hat{r}$$


But how about the force on object j? To do so, we invoke Newton's third law: Every force has an equal and opposite force:
$$(1.10): \vec{F_{i}} = -\vec{F_{j}}$$

Using equation 1.10, we have:
$$(1.11): \vec{F_{j}} = \frac{G M{i} M{j}}{r^{2}} \hat{r}$$

We can make the attractive nature of gravity more explicit by adding the subscripts i and j. By defining the vector relationship for any vector V:
$$(1.12): \vec{V_{ji}} = \vec{V_{i}} - \vec{V_{j}}$$

We can rewrite equation 1.2, the force of object j on object i, as:
$$(1.13): \vec{F_{ji}} = \frac{-G M{i} M{j}}{r_{ji}^2} \hat{r_{ji}}$$

And the reciprocal force, in equation 1.10 as:
$$(1.14): \vec{F_{ij}} = \frac{-G M{i} M{j}}{r_{ij}^2} \hat{r_{ij}}$$

 

[Sagekilla]

Solving the Equations of Motion

Where do we go from here? We can calculate, for any object, the acceleration due to the gravitational force of any other object. How can we calculate how a pair of particles move through space?

To do so, we need to make use of differential equations. To start with, we know some basic relationships: Velocity is the change in position with time, and acceleration is the change in velocity over time.

In mathematical form, we have:
$$(2.1): \vec{v} = \frac{d \vec{r}}{dt}$$
$$(2.2): \vec{a} = \frac{d \vec{v}}{dt} = \frac{d^{2} \vec{r}}{dt^{2}}$$

We can state Newton's second law as:
$$(2.3) \vec{F} = m \frac{d^{2} \vec{r}}{dt^{2}}$$

In the case of gravity, we obtain the following differential equation:
$$(2.4): M_{i} \frac{d^{2} \vec{r}}{dt^{2}} = \frac{-G M_{i} M_{j}}{r_{ji}^{2}}\hat{r_{ji}}$$
$$(2.5): \frac{d^{2} \vec{r}}{dt^{2}} = \frac{-G M_{j}}{r_{ji}^{2}}\hat{r_{ji}}$$

What does this tell us? It says that the second derivative of position is equal to some function of position. But how can we possibly hope to integrate this equation and find the position? Kepler was able to solve this problem analytically for two objects, but as we move to larger systems of objects it becomes impossible.

Our only solution now is to numerically find a solution to this differential equation.

Recall the definition of a derivative:
$$(2.6): \frac{dr}{dt} = \lim_{h \to 0} \frac{r(t + h) - r(t)}{h}$$

If we treat equation 2.5 as a difference equation, where we approximate the derivative by using a very small (but finite, and not infinitesimally small) stepsize dt, then we have:

$$(2.7): \frac{dr}{dt} \approx \frac{r(t + dt) - r(t)}{dt}$$

Rearranging terms in equation 2.6, we have:
$$(2.8): r(t) + \frac{dr}{dt} dt \approx r(t + dt)$$

Remember from equation 2.1 though, that velocity is the time derivative of position:
$$(2.8): r(t) + v(t) dt \approx r(t + dt)$$

We now have a way to solve for position at some new time given the previous position and velocity. If we apply the same kind of analysis to velocity, we obtain:
$$(2.9): v(t) + a(t) dt \approx v(t + dt)$$

If we take these to be true in the limiting case as dt approaches zero (Which will not be proven here -- You can find proofs online that explain this better than I can):
$$(2.10): r(t + dt) = r(t) + v(t) dt$$
$$(2.11): v(t + dt) = v(t) + a(t) dt$$

The method used to find equations 2.10 and 2.11 is known as the Euler method. If we perform some error analysis using a taylor series (Wikipedia has a wonderful article on this) we can see that the error will be proportional to the

square of dt, which puts it as a first order method. There are higher order methods which introduce less error per time step as well as special symplectic methods which preserve energy (Which the Euler method does not).

 

[Sagekilla]

Simulating Gravity in an N-body Simulation

So now we can solve the equations of motion for a system of bodies under mutual gravitational attraction. Our three equations to remember are:
$$(1.13): \vec{F_{ji}} = \frac{-G M{i} M{j}}{r_{ji}^2} \hat{r_{ji}}$$
$$(2.10): \vec{r}(t + dt) = \vec{r}(t) + \vec{v}(t) dt$$
$$(2.11): \vec{v}(t + dt) = \vec{v}(t) + \vec{a}(t) dt$$

To simulate this, we can follow this procedure:

1) Calculate the pairwise interactions between each body
2) Integrate the position of each body using equation 2.10
3) Integrate the velocity of each body using equation 2.11
4) Repeat for as many steps as you'd like.

Here's some sample C++ code that will perform the process above. Try to experiment a bit. Change how the particles are initialized, or create elaborate setups yourself.
Be aware the code below will not print anything out! If you want output, try opening a file and outputting the state of the whole system at every time step. I recommend www.cplusplus.com[/URL] as a reference.

[CODE]
#include <cmath>
#include <cstdlib>
#include <ctime>

using namespace std;

// A helper structure
struct Particle
{
float mass;
float pos[3];
float vel[3];
float acc[3];
};

void initParticles(int numParticles, Particle *p);
void calcGavity(int numParticles, Particle *p);
void integrate(int numParticles, Particle *p, float dt);

int main()
{
// Randomize starting seed
srand((unsigned) time(0));

// Set number of particles, step size, and steps
const int numParticles = 100;
const float dt = 0.01;
const int steps = 1000;

// Reserve space for numParticles
Particle *p = new Particle[numParticles];

// Initialize the particles
initParticles(numParticles, p);

for (int t = 0; t < steps; t++)
{
// Calculate pairwse interactions
// due to gravity
calcGravity(numParticles, p);

// Integrate the equations of motion
// using the Euler method (inaccurate!)
integrate(numParticles, p, dt);
}

// Free up memory reserved for the particles
delete[] p;

return 0;
}

void initParticles(int numParticles, Particle *p)
{
// Initialize each particle with a random
// position and velocity in range [-0.5, 0.5]
// and mass 1 / numParticles, for total mass = 1
for (int i = 0; i < numParticles; i++)
{
p[i].mass = 1.0 / numParticles;

for (int j = 0; j < 3; j++)
{
p[i].pos[j] = rand() / (float) RAND_MAX - 0.5;
p[i].vel[j] = rand() / (float) RAND_MAX - 0.5;
}
}
}

void calcGravity(int numParticles, Particle *p)
{
for (int i = 0; i < numParticles; i++)
{
// Clear out the acceleration for the ith particle
for (int k = 0; k < 3; k++)
p[i].acc[k] = 0;

for (int j = 0; j < numParticles; j+)
{
// Don't calculate gravity on yourself
if (i == j)
continue;

// Our position vector pointing from j to i
// And our distance between the two
float r[3] = {0, 0, 0};
float dist = 0;

// Get the vector from j to i
// and calculate distance squared
for (int k 0; k < 3; k++)
{
r[k] = p[i].pos[k] - p[j].pos[k];
dist += r[k] * r[k];
}

dist = sqrt(dist);
float distCube = dist * dist * dist;

// Calculate gravity: aji = -G * Mj * rji / |rji|^3
for (int k = 0; k < 3; k++)
p[i].acc[k] += -G * p[j].mass * r[k] / distCubed;
}
}
}

void integrate(int numParticles, Particle *p, float dt)
{
// Integrate using Euler's method:
// r(t + dt) = r(t) + v(t) * dt
// v(t + dt) = v(t) + a(t) * dt

for (int i = 0; i < numParticles; i++)
{
for (int j = 0; j < 3; j++)
{
p[i].pos[j] += p[i].vel[j] * dt;
p[i].vel[j] += p[i].acc[j] * dt;
}
}
}
[/CODE]

 

[Sagekilla]

Well, there's the tutorial! I hope this is helpful for all those who have a fascination with computers, physics, and simulations (as I do).

If anyone has any comments, suggestions, or criticism feel free to post it or PM me. I will make changes as necessary to correct any errors I may have made.