# Center of mass chaotic in n-body gravitational interactions?

I wrote a gravitational particle simulation. To prevent the particles flinging off out of the screen, I added a bit of friction (it also makes the animation prettier as particles spiral into one another :)). Total friction applied to particles is 0,05% per step or roughly 26% per second.

Since the particles start out scattered in random locations, I expected the center of mass to move, but not chaotically. Crashing of particles is accurately simulated. Can this behaviour seen in real life or is it a result of friction or some other factor I didn't think of? This is the result of particle simulation:

Green is the center of mass, white are individual particles.

EDIT:

To answer some of your questions, the particles can continue to go on forever out of the image. They will not be displayed, but they affect and are affected by every other particle.

My particles are attracted to each other by applying the following force to them at every step like this:

For $k$-th particle's $d$-th dimension at $n+1$-st step with mass $m$ where $\Delta l_d$ is the distance between particles in $d$-th dimension, $\Delta l$ is the total distance between particles, $D$ is the number of dimensions, and $\mu$ the factor of friction, velocity

$$v_{k,d,n+1}=(v_{k,d,n}+\sum_{i=1}^{k-1}{\frac{m_km_i\Delta l_d}{\Delta l^D}}+\sum_{i=k+1}^{n}{\frac{m_km_i\Delta l_d}{\Delta l^D}})(1-\mu)$$

I did some additional tests on the simulation. I removed friction, so there are no longer any "external" forces, and got this (note, the particle positions were generated using a random seed, so the starting conditions are not the same):

I bumped up the time step from 600 to 60000 per second, but it did not help. I figured it was rounding errors, so I switched from floating point numbers to much more accurate data structures with 1000 digits of precision. As expected, the simulation time became exponentially longer, but it did not really help. You can see sped up animation of high precision simulation here. I thought it may be the result of the fact that collisions have no elasticity, but thinking about it, it should not change the center of mass. Is it really the problem of still too small time step?

• I am no expert but n-body system with n above 2 are chaotic, so you expect multiple numerical errors to accumulate over time. In fact it is very common in these numerical simulations to correct the motions so that the CM stays always at the center of your simulation. – user83548 Nov 9 '15 at 14:05
• Yes, this is actually a hard problem to do in simulation correctly. What integration scheme are you using? – ja72 Nov 9 '15 at 14:08
• When a particle flies of the screen does it bounce back, wrap around or go on forever? – ja72 Nov 9 '15 at 14:09
• Presumably there is an exact equation, including the friction terms, that you are numerically modelling. If you could write that for us we might be able to see whether it conserves momentum. Two particles might be enough. – Keith McClary Nov 9 '15 at 22:52

According to conservation of momentum, the center of mass of a system cannot accelerate without external forces. In other words, if the center of mass starts out at rest (which is generally a good procedure in simulations), then it should always stay at rest.

It is normal for numerical errors to introduce deviations, but the motion you are seeing looks very large (really cool image by the way!). Introducing friction (effectively an external force) can also introduce deviations from momentum-conservation --- depending on how it is implemented. Both of these effects (friction and numerical error) can introduce chaotic perturbations to the center of mass.

The first thing you might want to try is seeing if decreasing your timestep makes the center of mass stay more stationary. The traditional criteria for determining what (maximum) timestep to use in hydrodynamic simulations is the "CFL Condition"; I have seem people generalize this to N-Body simulations by using the maximum particle velocity instead of the soundspeed, and the minimum separation between particles instead of the grid-size. This results in a condition like,

$$\Delta t \leq \frac{\Delta x_\textrm{min}}{v_\textrm{max}}$$

This condition can make the simulation become very 'expensive' (slow timesteps) when particles come close together - and their mutual forces (accelerations, and velocities) are very large. Often, to alleviate this, people will introduce 'gravitational softening', where the gravitational force between two-bodies is calculated as,

$$F = G \frac{m_1 m_2}{\left(r + a\right)^2}$$

where $a$ is the 'softening length', and should be small compared to the 'size of interest' (e.g. the initial separation of particles).

One quick and dirty fix for center of mass motion, is to calculate the velocity of the center of mass at each timestep, and subtract that velocity from all particles. This effectively re-centers the simulation on the center of mass at each timestep.

I assume that friction is an external velocity dependent force in your simulation code. Since you have such external forces, your total energy, total angular momentum, total momentum are likely not to be conserved.

In your case, the friction is a phenomenological external force, but similar behavior could also be simulated with a large particle, moving in a cloud of smaller particles. Gradually, the kinetic energy of this particle would be distributed equally among the degrees of freedom of this system. In this case, the total energy and the total momentum of the system will be conserved, since there are no external forces.

In general, I would recommend taking it step by step on numerical implementations, if correctness is required. Does single particle conserve energy? Does 'earth' orbiting around fixed 'sun' conserve energy and angular momentum? And if 'sun' is not fixed? When moving to chaotic systems, you could calculate the sum of total energy and the work done by the frictional force, and see if this is conserved?